In [44]:
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[44]:
In [2]:
%%html
<style>
td {
font-size: 14px
}
table {
float:left
}
</style>
NSBH Spin Analysis¶
In this work, I consider mass ratio 5 NSBH mergers with varying black hole spin. Neutron star is non-spinning with compactness of 0.15. Black Hole spin varies from -0.5 to 0.75. As a test study, the initial separation is kept small (only 7M). Below table provides additional simulation details.
| Simulation | Type | Mass Ratio | Separation | BH Spin | Resolution TOV | Resolution BH1 |
|---|---|---|---|---|---|---|
| Sim1 | NSBH | 5 | 7M | -0.5 | $R_*$/51 | $M_{BH}$/90 |
| Sim2 | NSBH | 5 | 7M | 0.0 | $R_*$/50 | $M_{BH}$/89 |
| Sim3 | NSBH | 5 | 7M | 0.25 | $R_*$/50 | $M_{BH}$/89 |
| Sim4 | NSBH | 5 | 7M | 0.5 | $R_*$/50 | $M_{BH}$/90 |
| Sim5 | NSBH | 5 | 7M | 0.75 | $R_*$/50 | $M_{BH}$/89 |
In [3]:
#Import Statements
import numpy as np
import matplotlib.pyplot as plt
import os, glob
from itertools import cycle
from shutil import copyfile
from scipy.interpolate import splrep, splev
import yt
from yt.visualization.eps_writer import multiplot_yt
from mpl_toolkits.axes_grid1 import AxesGrid
yt.funcs.mylog.setLevel(50)
import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)
In [4]:
from scipy.constants import G, c
Msun = 1.989e30
m1_in_Msun = c**2/G/Msun
s1_in_Msun = c**3/G/Msun
#print("1m = %g Msun, 1s = %g Msun" %(m1_in_Msun, s1_in_Msun))
In [5]:
#Make line plots, log and semilog plots of one variable
def lineplots(variable_name, filename, grid_coord='x',xlimit=None, ylimit=None, shiftx=0.,shifty=0., logx=False, logy=False, num_snapshots=None, title=None ):
''' Function - Make line plots, log and semilog plots of one variable
variable - string input - variable name
filename - string input - filepath
xlimit - optional - float - limit of x axis
logx - binary (True/False) - to change scale of x axis to log
logy - binary (True/False) - to change scale of y axis to log
title - title of the plot
num_snapshots - number of snapshots in time of variable
shiftx - shift x coordinate by this amount
'''
if grid_coord=='x': grid_coln = 9
elif grid_coord=='y': grid_coln = 10
else: grid_coln=11
t, x, var = np.loadtxt(filename, unpack=True, usecols=(8,grid_coln,12))
t_uniq = np.unique(t)
plt.figure(figsize=(15,8))
if num_snapshots==None:
for time_instance in t_uniq:
xx_instance = x[t==time_instance]
var_instance = var[t==time_instance]
idx_sorted = np.argsort(xx_instance)
xx_sort, var_sort = xx_instance[idx_sorted], var_instance[idx_sorted]
legend=False
plt.plot(xx_sort,var_sort)
else:
legend=True
time_snapshot_idx = np.linspace(0,len(t_uniq)-2, num_snapshots)
time_snapshot_idx = time_snapshot_idx.astype(int)
time_snapshot = t_uniq[time_snapshot_idx]
for time_instance in time_snapshot:
idx_instance = np.where(t==time_instance)
xx_instance = x[t==time_instance]
var_instance = var[t==time_instance]
idx_sorted = np.argsort(xx_instance)
xx_sort, var_sort = xx_instance[idx_sorted], var_instance[idx_sorted]
plt.plot(xx_sort,var_sort, label="t=%g"%time_instance)
plt.xlabel('x')
plt.ylabel(variable_name)
if not(title==None):
plt.title(title, fontsize=16)
if isinstance(xlimit,float):
plt.xlim(-1*xlimit+shiftx,xlimit+shiftx )
if isinstance(ylimit,float):
plt.ylim(-1*ylimit+shifty,ylimit+shifty )
if(logx==True): plt.xscale('log')
if(logy==True): plt.yscale('log')
if legend==True:
plt.legend()
plt.show()
plt.close()
def linesubplots(numvars, filename, yaxisname, labelname):
''' Function - Make subplots of multiple variables
numvars - integer - number of variables
filename - list of file paths
yaxisname - list of y-axis labels
labelname - list of labels
'''
num_rows, num_cols = 0, 0
if numvars%2==0:
num_rows, num_cols = numvars/2, 2
else:
num_rows, num_cols = numvars, 1
fig = plt.figure(figsize=(12,7*num_rows))
time, x, variable = {}, {}, {}
#Read data and store in dictionaries
for i in range(0,numvars):
tt, xx, data = np.loadtxt(filename[i], unpack=True, usecols=(8,9,12))
idx_sorted = np.argsort(xx)
xx, data = xx[idx_sorted], data[idx_sorted]
time['var%d'%i] = tt
x['var%d'%i] = xx
variable['var%d'%i] = data
plt.subplot(num_rows, num_cols, i+1)
plt.plot(x['var%d'%i], variable['var%d'%i], label=labelname[i])
plt.ylabel(yaxisname[i])
plt.xlabel('x')
plt.xlim(-8,8)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
def comparisonplots(varname, numdatasets, filename, labelname, xlimit=None, logx=False, logy=False, freq=1):
''' Function - single plots showing comparison of two or more variables
varname - string - variable name which is being compared
numdatasets - integer - number of datasets
filename - list of file paths
labelname - list of label of each dataset
'''
time, x, variable = {}, {}, {}
lines = ["-","--","-.",":"]
linecycler = cycle(lines)
fig = plt.figure(figsize=(15,8))
#Read data and store in dictionaries
for i in range(numdatasets):
tt, xx, data = np.loadtxt(filename[i], unpack=True, usecols=(8,9,12))
t_arr = np.unique(tt)
linestyle=next(linecycler)
for t_elem in t_arr:
xx_t = xx[np.where(tt==t_elem)]
data_t = data[np.where(tt==t_elem)]
idx_sorted = np.argsort(xx_t)
xx_t, data_t = xx_t[idx_sorted], data_t[idx_sorted]
plt.plot(xx_t, data_t, ls=linestyle, label=labelname[i])
plt.ylabel(varname)
plt.xlabel('x')
if(logx==True): plt.xscale('log')
if(logy==True): plt.yscale('log')
if isinstance(xlimit,float):
plt.xlim(-1*xlimit,xlimit)
plt.legend()
plt.show()
plt.close()
def comparisonmultiplots(varnames, filenames, labelnames, xlimit=None, logx=False, logy=False):
''' Create multiple subplots comparing two quantities in each subplot
varnames: list of name of each variable (this name is also used to label y axis)
numdatasets: number of datasets for each variable
filename: dictionary of filenames for each variable and each dataset. Use index in form of rows and columns where
each row contains all datasets for one subplot
labelname: list of names of labels for each varname. Should have size of numdatasets
'''
num_vars = len(varnames)
num_datasets = len(labelnames)
num_files = len(filenames.keys())
assert (num_files==num_vars*num_datasets), "Total number of filenames do not match with number of variables and labelnames"
if xlimit==None:
xlimit=30
num_cols=2
if num_vars%2==0:
num_rows = num_vars/2
else:
num_rows, num_cols = num_vars, 1
time, x, variable={},{},{}
plt.figure(figsize=(15, 8*num_rows))
linestyles = ["-","--","-.",":"]
for i in range(num_files):
tt, xx, data = np.loadtxt(filenames[sorted(filenames.keys())[i]], unpack=True, usecols=(8,9,12))
idx_sorted = np.argsort(xx)
xx, data = xx[idx_sorted], data[idx_sorted]
time[i] = tt
x[i] = xx
variable[i] = data
if i%num_datasets==0:
plt.subplot(num_rows, num_cols, 1+i/num_datasets)
plt.ylabel(varnames[i/num_datasets])
plt.xlabel('x')
plt.xlim(-1*xlimit, xlimit)
#print i
plt.plot(x[i], variable[i], ls = linestyles[i%num_datasets], label=labelnames[i%num_datasets])
plt.legend(fontsize=14)
plt.show()
plt.close()
In [6]:
os.chdir('/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_5/SpinStudy/')
In [7]:
#Directories
nsbh_q5_d7_a1_n0d5_dir1 = 'NSBH_q5_D7_a1_-0.5_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_0d0_dir2 = 'NSBH_q5_D7_a1_0.0_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_p0d25_dir3 = 'NSBH_q5_D7_a1_0.25_a2_0.0_C0.15_M107/'
nsbh_q5_d7_a1_p0d5_dir4 = 'NSBH_q5_D7_a1_0.5_a2_0.0_C0.15_M109/'
nsbh_q5_d7_a1_p0d75_dir5 = 'NSBH_q5_D7_a1_0.75_a2_0.0_C0.15_M107/'
Madm_nsbh_a1_0d0 = 0.9938816335979512306
Madm_nsbh_a1_p0d25 = 0.9937419751277200008
Madm_nsbh_a1_p0d5 = 0.993621911714029471
Madm_nsbh_a1_p0d75 = 0.9935297248578904838
Madm_nsbh_a1_n0d5 = 0.9942234586460986234
#Initial Configuration
Dist_from_src=100. #Mpc
Mtotal = 1.35*3 #Msun
mass_ratio=5
MPc = 3.086e+22 #metres
Mtotal_to_MPc = (Mtotal/m1_in_Msun/MPc)
nsbh_a1_n0d5 = {'dir':nsbh_q5_d7_a1_n0d5_dir1, 'label':r'$a_{BH} = -0.5$', 'type':'NSBH', 'symmetry':True, \
'ls':'-', 'color':'k', 'M_adm':Madm_nsbh_a1_n0d5} #, 'tidal_radius':6.6}
nsbh_a1_0d0 = {'dir':nsbh_q5_d7_a1_0d0_dir2 , 'label':r'$a_{BH} = 0$', 'type':'NSBH', 'symmetry':True, \
'ls':':', 'color':'b', 'M_adm':Madm_nsbh_a1_0d0}
nsbh_a1_0d25 = {'dir':nsbh_q5_d7_a1_p0d25_dir3, 'label':r'$a_{BH} = 0.25$' ,'type':'NSBH', 'symmetry':True, \
'ls':':', 'color':'orangered', 'M_adm':Madm_nsbh_a1_p0d25}
nsbh_a1_0d5 = {'dir':nsbh_q5_d7_a1_p0d5_dir4, 'label':r'$a_{BH} = 0.5$', 'type':'NSBH', 'symmetry':True, \
'ls':':', 'color':'m', 'M_adm':Madm_nsbh_a1_p0d5}
nsbh_a1_0d75 = {'dir':nsbh_q5_d7_a1_p0d75_dir5, 'label':r'$a_{BH} = 0.75$', 'type':'NSBH', 'symmetry':True, \
'ls':':', 'color':'forestgreen', 'M_adm':Madm_nsbh_a1_p0d75}
compare_nsbh = {'nsbh1':nsbh_a1_n0d5, 'nsbh2':nsbh_a1_0d0, 'nsbh3':nsbh_a1_0d25, 'nsbh4':nsbh_a1_0d5, 'nsbh5':nsbh_a1_0d75}
figdir = '/localdata2/bkhamesra3/simulations/Hive/NSBH/Spin_Comparison/Sep_7M'
if not os.path.exists(figdir):
os.makedirs(figdir)
Merger Time¶
- $T_{Peak}$ - Time of $\psi_4$ peak amplitude.
- $T_{Peak} - R_{ex}$ - Approximate Retarded Time of peak amplitude defined as $t' = t - R_{ex}$. This is used as the definition of the Time of Merger in this analysis.
- $T_{Peak}$ (Retarded) - Retarded Time of peak amplitude defined as - $t' = t - R_{ex} - 2*log(\frac{R_{ex}}{2}-1)$
In [8]:
# Define Merger Time in retarded frame - t' = t - D - 2*log(D/2-1) -
# Reference - Equation 26, https://iopscience.iop.org/article/10.1088/0264-9381/31/2/025012/pdf
print("Simulation \t \tT_peak \t T_peak - R_ex \t T_peak(retarded)")
print("-"*100)
for i, key in enumerate(sorted(compare_nsbh.keys())):
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
psi4 = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4):
psi4 = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t, re, im = np.loadtxt(psi4, unpack=True, usecols=(0,1,2))
det_radius = float(((psi4.split('_')[-1]).split('.asc')[0]).split('r')[1])
amp = np.sqrt(re**2 + im**2)
tmax_det = t[amp==np.max(amp)]
tmax_ret = tmax_det-det_radius -2.*np.log(det_radius/2. - 1.)
tmerger_src = tmax_det-det_radius
compare_nsbh[key]['t_retarded'] = tmax_ret
compare_nsbh[key]['t_merger_src'] = tmerger_src
print("%15s \t %.2f \t %.2f \t %.2f" %(compare_nsbh[key]['label'],tmax_det, compare_nsbh[key]['t_merger_src'],compare_nsbh[key]['t_retarded']))
In [9]:
#Compute ADM angular momentum
for i, key in enumerate(sorted(compare_nsbh.keys())):
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
parameter_file = glob.glob(os.path.join(datadir, '*.par'))[0]
labelname = compare_nsbh[key]['label']
#adm_angmom = compare_nsbh[key]['J_adm']
x1,y1,z1,x2,y2,z2,px1,py1,pz1,px2,py2,pz2,jx1,jy1,jz1,jx2,jy2,jz2 = np.zeros(18)
with open(parameter_file) as parfile:
for line in parfile.readlines():
if 'BowenID::object_rx[0]' in line:
x1 = float(line.split('=')[1])
elif 'BowenID::object_rx[1]' in line:
x2 = float(line.split('=')[1])
elif 'BowenID::object_Px[0]' in line:
px1 = float(line.split('=')[1])
elif 'BowenID::object_Py[0]' in line:
py1 = float(line.split('=')[1])
elif 'BowenID::object_Pz[0]' in line:
pz1 = float(line.split('=')[1])
elif 'BowenID::object_Px[1]' in line:
px2 = float(line.split('=')[1])
elif 'BowenID::object_Py[1]' in line:
py2 = float(line.split('=')[1])
elif 'BowenID::object_Pz[1]' in line:
pz2 = float(line.split('=')[1])
elif 'BowenID::object_Jx[0]' in line:
jx1 = float(line.split('=')[1])
elif 'BowenID::object_Jy[0]' in line:
jy1 = float(line.split('=')[1])
elif 'BowenID::object_Jz[0]' in line:
jz1 = float(line.split('=')[1])
elif 'BowenID::object_Jx[1]' in line:
jx2 = float(line.split('=')[1])
elif 'BowenID::object_Jy[1]' in line:
jy2 = float(line.split('=')[1])
elif 'BowenID::object_Jz[1]' in line:
jz2 = float(line.split('=')[1])
r_bh = np.array((x1, y1, z1))
p_bh = np.array((px1, py1, pz1))
j_bh = np.array((jx1, jy1, jz1))
r_ns = np.array((x2, y2, z2))
p_ns = np.array((px2, py2, pz2))
j_ns = np.array((jx2, jy2, jz2))
j_adm = np.cross(r_bh, p_bh ) + np.cross(r_ns, p_ns) + j_bh + j_ns
compare_nsbh[key]['J_adm'] = j_adm
Density Snapshots¶
We look at the density evolution of the NSBH systems. Figures from top to bottom are arranged in the increasing order of the BH spins. For each case, the panels on the top from left to right show time snapshots at 100M, 50M and 20M before the merger while panels on the bottom show the snapshots at the merger, and 100M and 300 M after the merger. A clear impact of the merger is the increasing disruption of the star which is minimal for $a_{BH} = -0.5$ and maximum for $a_{BH} = 0.75$.
In [10]:
#Define Fields for YT
from yt import derived_field
def _density_gu(field, data):
'''Specifies density in geometrical units'''
return data["density"]
def _density_cgs_q5(field, data):
'''Specifies density in geometrical units'''
Mt = 1.35*6
return (data["density"]/Mt**2)
#yt.add_field(("gas", "density_gu"), function=_density_gu, units="g/cm**3")
yt.add_field(("gas", "density_gu"), function=_density_gu, units="code_mass/code_length**3")
#yt.add_field(("gas", "density_cgs_q2"), function=_density_cgs_q2, units="g/cm**3")
#yt.add_field(("gas", "density_cgs_q3"), function=_density_cgs_q3, units="g/cm**3")
yt.add_field(("gas", "density_cgs_q5"), function=_density_cgs_q5, units="g/cm**3")
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
#print(labelname+'\n')
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merger_src = compare_nsbh[key]['t_merger_src']
# print("%s:Time of merger = %g"%(labelname, t_merger_src))
else:
continue
M_in_Msun = 1.35*6
fieldname = "density_cgs_q5"
#filename = glob.glob(os.path.join(datadir,'../../output-0000/NSBH*/rho.x.asc'))[0]
#t, x, rho = np.loadtxt(filename, unpack=True, usecols=(8,9,12))
#x = x[t==0]
#rho = rho[t==0]
max_rho = 0.101473916805118 #np.amax(rho)
#x_maxrho = x[np.argmax(rho)]
max_rho_si = (max_rho/(M_in_Msun**2))*(Msun*m1_in_Msun**3)
h5file = os.path.join(datadir, 'rho.xy.h5')
varname1 = fieldname
ds1 = yt.load(h5file, iteration=0) # load data
_, max_rho_si_yt = ds1.all_data().quantities.extrema(varname1)
varname2 = 'density_gu'
ds2 = yt.load(h5file, iteration=0) # load data
min_rho_gu_yt, max_rho_gu_yt = ds2.all_data().quantities.extrema(varname2)
#print("t=0: x=%g, rho_c = %g/M^2 (%g g/cm^3), rho_c (yt) = %g /M^2 (%g g/cm^3)"%(x_maxrho, max_rho, max_rho_si/1e3, max_rho_gu_yt, max_rho_si_yt))
In [11]:
#Rest Mass Density Contour Plots
import yt
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
mpl.rc('axes', labelsize=16, grid=False)
# See http://matplotlib.org/mpl_toolkits/axes_grid/api/axes_grid_api.html
# These choices of keyword arguments produce a four panel plot that includes
# four narrow colorbars, one for each plot. Axes labels are only drawn on the
# bottom left hand plot to avoid repeating information and make the plot less
# cluttered.
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
spin = float(labelname.split('=')[1][:-1])
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merger_src = compare_nsbh[key]['t_merger_src']
else:
continue
fig = plt.figure(figsize=(10,8))
#Setup axes grid
grid = AxesGrid(fig, (0.1,0.1,1.2,1.2),
nrows_ncols = (2, 3),
axes_pad = 0.5,
label_mode = "L",
share_all = False,
cbar_location="right",
cbar_mode="single",
cbar_size="3%",
cbar_pad="0%")
#Define variable name for plots
varname = 'density'
#Path for density and BH diagnostics data
h5file = os.path.join(datadir, 'rho.xy.h5')
bhdiag = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
shifttracker = os.path.join(datadir, 'ShiftTracker0.asc')
assert(os.path.exists(h5file)),"H5 file not found"
M = 1.35*6
#varname_cgs = "density_cgs_q5"
#Read data from BH diagnostics file
it_st, t_st = np.loadtxt(shifttracker, unpack=True, usecols=(0,1))
it, t, xbh, ybh, zbh, rbh, bh_mass = np.loadtxt(bhdiag, unpack=True, usecols=(0,1,2,3,4,7,26))
idx_merger = np.amin(np.where(t>=t_merger_src))
#print(t_merger_src, t[idx_merger])
t_centered = t - t[idx_merger]
t_st_centered = t_st - t[idx_merger]
t_premerger1 = -100
t_premerger2=-50
t_premerger3=-20
t_postmerger1=0
t_postmerger2=100
t_postmerger3=300
time_inst = [t_premerger1, t_premerger2, t_premerger3, t_postmerger1, t_postmerger2, t_postmerger3]
ds_allit = yt.load(h5file)
iterations = ds_allit.get_all_iterations()
k=0
print("BH Spin = {} \n\t Pre Merger - t1 = {}({}), t2 = {}({}), t3 = {}({}),\
\n\t Post Merger - t1 = {}({}), t2 = {}({}), t3 = {}({})\n"\
.format(spin, t_premerger1,t_premerger1 + t[idx_merger], t_premerger2, t_premerger2 + t[idx_merger],\
t_premerger3, t_premerger3 + t[idx_merger], t_postmerger1, t_postmerger1 + t[idx_merger],t_postmerger2,\
t_postmerger2 + t[idx_merger],t_postmerger3, t_postmerger3 + t[idx_merger]))
for t_pm in time_inst:
it_pm_st = it_st[np.amax(np.where(t_st_centered<=t_pm))]
idx_pm = np.amax(np.where(t_centered<=t_pm))
it_pm = it [idx_pm]
xbh_pm = xbh[idx_pm]
ybh_pm = ybh[idx_pm]
zbh_pm = zbh[idx_pm]
rbh_pm = rbh[idx_pm]
iter_premerger = np.amax(np.where(iterations<=it_pm))
ds = yt.load(h5file, iteration=iterations[iter_premerger]) # load data
min_var, max_var = ds.all_data().quantities.extrema(varname)
#As we cannot change units while loading the data, we call for new variables with specific units
#min_rho_gu, max_rho_gu = ds.all_data().quantities.extrema('density_gu')
#min_rho_si, max_rho_si = ds.all_data().quantities.extrema(varname_cgs)
#print("t=%d: Maximum Density = %g/M^2 (%g g/cm^3)"%(t_pm, max_rho_gu, max_rho_si))
axis_width=30
#Plotting density is fine here as we later set the units to code units which plots exact data stored in the file
p = yt.SlicePlot(ds, 'z', varname, axes_unit='code_length', origin='native', width=axis_width, center = [0,0,0])
p.set_unit('density', 'code_mass/code_length**3')
p.set_font_size(16)
p.set_xlabel('x/M')
p.set_ylabel('y/M')
p.set_colorbar_label('density','Density (1/M$^2$)')
p.set_cmap(varname, 'dusk') #'viridis')
p.annotate_sphere([xbh_pm, ybh_pm, zbh_pm], radius=(rbh_pm, 'code_length'), circle_args={'facecolor':'black','fill':1, 'edgecolor':'white'})
#p.annotate_contour(varname, plot_args={'colors':'white', 'linestyles':'dotted'})
p.set_font_size(20)
#p.annotate_text((0.45,0.9), 't = %.2f'%t[idx_pm], coord_system='figure', text_args={'color':'white'})
#if k==i:p.annotate_title('q=%d'%mass_ratio)
p.set_zlim(varname, 0.05, 1e-7)
plot = p.plots[varname]
plot.figure = fig
plot.axes = grid[k].axes
plot.cax = grid.cbar_axes[k]
# Actually redraws the plot.
p._setup_plots()
k+=1
#plt.savefig( os.path.join(figdir, 'DensitySnapshots_q%d.png')%mass_ratio, dpi=100 )
plt.show()
Density Oscillations__¶
Below figure shows the density oscillations in the neutron star and corresponding frequency of oscillations. The black hole spin barely impacts the NS oscillations in amplitude or frequency.
In [12]:
#Density Oscillations and FT Plots - Thesis figures
i=0
filename = 'MinSearch0.asc'
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(15,7))
#color = [ 'k','b','orangered','k','b','orangered' ,'k','k','k','brown','m','coral','forestgreen']
#ls = ['-', '--', ':','-','-','-','--',':','-.']
for i, direc in enumerate(sorted(compare_nsbh.keys())):
if compare_nsbh[direc]['type']=='BBH':
continue
labelname = compare_nsbh[direc]['label']
datadir = os.path.join(compare_nsbh[direc]['dir'], 'Summary/Data')
color=compare_nsbh[direc]['color']
filepath = os.path.join(datadir, filename)
t, max_rho = np.loadtxt(filepath, unpack=True, usecols=(1,5))
idx_stablestar = np.intersect1d(np.where(t>50), np.where(t<150))
mtotal = 1.35*(1+5)
idx_fft = np.intersect1d(np.where(t>0), np.where(t<150))
rho_fftfreq = np.fft.fftfreq(np.size(max_rho[idx_fft]), d=(t[2]-t[1]))
rho_fft = np.fft.fft(max_rho[idx_fft]/max_rho[0])
ax1.plot(t, (max_rho)/max_rho[0], ls='-',c=color, label=labelname)
ax1.set_xlabel('T/M')
ax1.set_ylabel(r'$ \rho_{max}/\rho_{init}$')
ax1.legend()
#ax1.set_ylim(0.88, 1.5)
ax1.set_xlim(0,150)
ax1.grid(True)
ax2.plot(rho_fftfreq/(mtotal/s1_in_Msun), rho_fft.real, ls='-',c=color, label=labelname)
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel(r'$FFT[\rho_{max}/\rho_{init}]$')
ax2.set_ylim(-50,70)
ax2.set_xlim(500,2000)
ax2.legend()
ax2.grid(True)
i+=1
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'ReviewFigures/RevFig_DensityOscillations.png'), dpi=100)
plt.show()
plt.close()
Orbital Dynamics¶
Below is comparison of trajectories and orbital separation. Note that the initial momenta values differ between each case which can lead to differences in eccentricity, trajectories and number of orbits. Few key observations -
- The merger time and number of orbits increase with spin and/or eccentricity.
- Higher spin leads to larger eccentricity as reflected in oscillations in the orbital separation and require better eccentricity reduction techniques.
- The BH velocity and COM offset also increases with spin of the BH at merger.
- System facing strongest disruption ($a_{BH} = 0.75$) have highest BH speed at the merger and also suffers maximum changes in the speed. As such, its COM offset remains roughly same even after merger.
In [13]:
#Binary trajectory parameters
mpl.rc('axes', labelsize=16, grid=True)
plt.figure(figsize=(15,20))
for i, key in enumerate(sorted(compare_nsbh.keys())):
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
label_obj2='NS'
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
label_obj2='BH'
t_postmerger = compare_nsbh[key]['t_merger_src']
label = compare_nsbh[key]['label']
bh1_search= 'ShiftTracker0.asc'
filepath = os.path.join(datadir, bh1_search)
if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
t1,x1,y1,z1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4))
ns2_search= 'volume_integrals-GRMHD.asc'
bh2_search= 'ShiftTracker1.asc'
if (compare_nsbh[key]['type']=='NSBH'):
filepath = os.path.join(datadir, ns2_search)
t2,x2,y2,z2,rest_mass = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4))
x2 = x2/rest_mass
y2 = y2/rest_mass
z2 = z2/rest_mass
elif(compare_nsbh[key]['type']=='BBH'):
filepath = os.path.join(datadir, bh2_search)
t2,x2,y2,z2 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4))
plt.subplot(3,2,i+1)
plt.plot(x1[t1<t_postmerger],y1[t1<t_postmerger], label='BH')
plt.plot(x2[t2<t_postmerger],y2[t2<t_postmerger], label=label_obj2)
plt.xlabel('x/M')
plt.ylabel('y/M')
plt.title(label, fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
Final BH Offset and Coordinate velocity at merger¶
- BH offset increases with spin
- Speed increases with spin
In [14]:
#Orbital Parameters
ls=['-','--',':','-.', '--',':']
c = ['r','b','k','m','coral','forestgreen']
orbdyn_dict = {}
print('-'*110)
print("System \t \t (T_peak - R)\t \t(x, y) of BH \t \t BH Offset \t (vx, vy) of BH \t Speed")
print('-'*110)
for i, key in enumerate(sorted(compare_nsbh.keys())):
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
label = compare_nsbh[key]['label']
t_merger = compare_nsbh[key]['t_merger_src']
bhsearch= 'ShiftTracker0.asc'
filepath = os.path.join(datadir, bhsearch)
if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
t1,x1,y1,z1,vx1,vy1,vz1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
x1 = x1 [t1<t_merger]
y1 = y1 [t1<t_merger]
z1 = z1 [t1<t_merger]
vx1 = vx1[t1<t_merger]
vy1 = vy1[t1<t_merger]
vz1 = vz1[t1<t_merger]
t1 = t1[t1<t_merger]
t_com, x_com, y_com, z_com, vx_com, vy_com, vz_com = t1[-1], x1[-1], y1[-1], z1[-1],\
vx1[-1], vy1[-1], vz1[-1]
com_offset = (x_com**2 + y_com**2 + z_com**2)**0.5
com_velocity = np.sqrt(vx_com**2 + vy_com**2 + vz_com**2)
print("%15s \t %.3g \t \t (%.2g, %.2g) \t%.3g \t \t(%.2g, %.2g) \t %.2g"%(label, t_com, x_com, y_com,\
com_offset,vx_com, vy_com, com_velocity))
ns2_search= 'volume_integrals-GRMHD.asc'
if (not os.path.exists(os.path.join(datadir, ns2_search))) and compare_nsbh[key]['type']=='NSBH':
raise NameError("%s:Volume Integrals file not found. Minsearch will be used. "%compare_nsbh[key]['label'])
bh2_search= 'ShiftTracker1.asc'
if (compare_nsbh[key]['type']=='NSBH'):
filepath = os.path.join(datadir, ns2_search)
t2,x2,y2,z2,rest_mass = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4))
x2 = x2[t2<t_merger]/rest_mass[t2<t_merger]
y2 = y2[t2<t_merger]/rest_mass[t2<t_merger]
z2 = z2[t2<t_merger]/rest_mass[t2<t_merger]
t2 = t2[t2<t_merger]
vx2 = np.gradient(x2)/np.gradient(t2)
vy2 = np.gradient(y2)/np.gradient(t2)
vz2 = np.gradient(z2)/np.gradient(t2)
elif (compare_nsbh[key]['type']=='BBH'):
filepath = os.path.join(datadir, bh2_search)
t2,x2,y2,z2,vx2,vy2,vz2 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
x2 = x2 [t2<t_merger]
y2 = y2 [t2<t_merger]
z2 = z2 [t2<t_merger]
vx2 = vx2[t2<t_merger]
vy2 = vy2[t2<t_merger]
vz2 = vz2[t2<t_merger]
t2 = t2[t2<t_merger]
t_common = np.intersect1d(t1, t2)
idx_st = np.isin(t1,t_common)
idx_ms = np.isin(t2,t_common)
sep = np.sqrt((x1[idx_st]-x2[idx_ms])**2 + (-y2[idx_ms]+y1[idx_st])**2 + (-z2[idx_ms]+z1[idx_st])**2)
orbphi = np.unwrap(np.arctan2(-y2[idx_ms]+y1[idx_st],-x2[idx_ms]+x1[idx_st]))
orbtheta = np.arccos(np.divide(z1[idx_st]-z2[idx_ms],sep))
vx = vx1[idx_st]-vx2[idx_ms]
vy = vy1[idx_st]-vy2[idx_ms]
vz = vz1[idx_st]-vz2[idx_ms]
#Velocities computed for BH are much more accurate as it does not use FD and computed from Shift. As we lack
#data for every iteration in NS + use second order FD, velocities tend to have higher numerical errors. Hence
#compute rdot and phidot directly from numerical derivative of the functions
#rdot = (vx*np.cos(orbphi) + vy*np.sin(orbphi))*np.sin(orbtheta) + vz*np.cos(orbtheta)
#phidot = np.divide((vy*np.cos(orbphi) - vx*np.sin(orbphi)), (sep*np.sin(orbtheta)))
rdot = np.gradient(sep)/np.gradient(t1[idx_st])
phidot = np.gradient(orbphi)/np.gradient(t1[idx_st])
orbdyn_dict[key] = {}
orbdyn_dict[key]['time'] = t_common
orbdyn_dict[key]['separation'] = sep
orbdyn_dict[key]['phase'] = orbphi
orbdyn_dict[key]['rdot'] = rdot
orbdyn_dict[key]['omega'] = phidot
Final BH Offset and Coordinate velocity at 300M after merger¶
- BH Offset continues to increase but the differences with offset at merger decrease with spin
- Coordinate speed drastically decreases for large spin values and do not follow the same trend as before.
In [15]:
#Orbital Parameters
ls=['-','--',':','-.', '--',':']
c = ['r','b','k','m','coral','forestgreen']
sep_list = []
orbphase_list = []
rdot_list = []
phidot_list = []
t_common_list = []
label_list = []
print('-'*110)
print("System \t \t (T_peak - R)+300 \t (x, y) of BH \t \tBH Offset \t (vx, vy) of BH \t Speed")
print("-"*110)
for i, key in enumerate(sorted(compare_nsbh.keys())):
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_postmerger = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_postmerger = compare_nsbh[key]['t_merger_src']
label = compare_nsbh[key]['label']
bhsearch= 'ShiftTracker0.asc'
filepath = os.path.join(datadir, bhsearch)
if not os.path.exists(filepath): raise RuntimeError("%s file not found"%filepath)
t1,x1,y1,z1,vx1,vy1,vz1 = np.loadtxt(filepath, unpack=True, usecols=(1,2,3,4,5,6,7))
if t1[-1]>t_postmerger:
idx_300 = np.amin(np.where(t1>=t_postmerger+300))
t_com, x_com,y_com, z_com, vx_com, vy_com, vz_com = t1[idx_300], x1[idx_300], y1[idx_300], z1[idx_300],\
vx1[idx_300], vy1[idx_300], vz1[idx_300]
com_offset = (x_com**2 + y_com**2 + z_com**2)**0.5
com_velocity = np.sqrt(vx_com**2 + vy_com**2 + vz_com**2)
if y_com<0:
print("%15s \t %.4g \t \t (%.2f, %.2f) \t%.3g \t\t (%.2g, %.2g) \t %.2g"%(compare_nsbh[key]['label'],\
t_com, x_com, y_com, com_offset, \
vx_com, vy_com, com_velocity ))
else:
print("%15s \t %.4g \t \t (%.2f, %.2f) \t%.3g \t\t (%.2g, %.2g) \t %.2g"%(compare_nsbh[key]['label'],\
t_com, x_com, y_com, com_offset, \
vx_com, vy_com, com_velocity ))
Orbital Separation and Phase¶
In [16]:
#Orbital Separation and Phase - NSBH vs BBH
fig, ax = plt.subplots(2,1,figsize=(10,10))
print('-'*100)
print("Simulation \t Merger Time \t Number of Orbits")
print('-'*100)
keys = compare_nsbh.keys()
for i,key in enumerate(sorted(compare_nsbh.keys())):
label = compare_nsbh[key]['label']
color = compare_nsbh[key]['color']
ls = compare_nsbh[key]['ls']
print("{:<16s} \t {:>4.3g} \t \t{:>4.3g} ".format(label,orbdyn_dict[key]['time'][-1], orbdyn_dict[key]['phase'][-1]/2./np.pi))
#Derivatives of mixed binary can be extremely noise. So instead of using the direct data, we perform
#fitting using spline interpolation and use them to compute the derivatives
sep_splinefit = splrep(orbdyn_dict[key]['time'],orbdyn_dict[key]['separation'],k=5,s=3)
phase_splinefit = splrep(orbdyn_dict[key]['time'], orbdyn_dict[key]['phase'],k=5,s=3)
orbdyn_dict[key]['sepfit'] = sep_splinefit
orbdyn_dict[key]['phasefit'] = phase_splinefit
if 'nsbh' in key:
#td_radius = compare_nsbh[key]['tidal_radius']
#Compute time of reaching tidal radius
time_key = orbdyn_dict[key]['time']
sep_key = orbdyn_dict[key]['separation']
#time_tidalrad = time_key[np.amin(np.argwhere(sep_key<=td_radius))]
ax[0].plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['separation'], ls='-',c=color,label=label)
#plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], sep_splinefit), ls='--',c='k', lw=1,label='Spline Fit')
ax[0].set_xlabel('T/M', fontsize=16)
ax[0].set_ylabel('Separation', fontsize=16)
ax[0].legend()
ax[1].plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['phase'], ls='-',c=color,label=label)
#plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], phase_splinefit), ls=':',c='k',label='Spline Fit')
ax[1].set_xlabel('T/M', fontsize=16)
ax[1].set_ylabel('Orbital Phase', fontsize=16)
ax[1].legend()
print('\n\n')
plt.tight_layout()
fig.suptitle('Figure 2.1 - NSBH vs BBH', fontsize=20, y=1.02)
plt.show()
plt.close()
Psi4 - (2,2) Mode Comparison¶
In [17]:
#Re(Psi4) plots:
import matplotlib.ticker as plticker
fig, ax = plt.subplots(1, 1, figsize=(15,7))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
color=compare_nsbh[key]['color']
ls = compare_nsbh[key]['ls']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
r = float(((filepath.split('.asc')[0]).split('_')[-1]).split('r')[-1])
ax.plot(t22, r*psi4_re22,c=color,ls='-', label=labelname)
#ax.plot(t22,amp,ls=ls,c=color)
ax.set_xlabel(r'$T/M$')
ax.set_ylabel(r'$Re(\psi_4)\,r\,M$')
ax.set_xlim(0,900)
ax.legend(loc=3)
ax.grid(True)
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax.xaxis.set_major_locator(loc)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'C3_Psi4Plots.png'), dpi=200)
plt.show()
plt.close()
Psi4 Amplitude¶
In [18]:
#Psi4 Amplitude plots: BBH vs NSBH
i=0
color = ['k','b', 'orangered','r','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
fig=plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
col = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls = '-'
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls=':'
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
R = float(((filepath.split('.asc')[0]).split('_')[-1]).split('r')[1])
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
psi4_re22 = psi4_re22[t22>100]
psi4_im22 = psi4_im22[t22>100]
t22 = t22[t22>100]
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t22-t22[max_amp_idx], R*amp,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T - R_{ex}(M)$')
plt.ylabel(r'$r \, |\Psi_4|M$')
plt.xlim(-50,100)
plt.ylim(1e-5, 0.03)
plt.legend()
plt.grid(True)
if "-0.5" in labelname:
continue
plt.subplot(122)
plt.plot(t22, R*amp,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$r \, |\Psi_4|M$')
plt.xlim(150,300)
plt.ylim(0,0.002)
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.grid(True)
plt.legend()
plt.tight_layout()
#fig.suptitle('Figure 2: BBH vs NSBH Psi4 (2,2) Mode', fontsize=16, y=1.08)
#plt.savefig(os.path.join(figdir,'C3_Psi4Amp.png'), dpi=100)
plt.show()
plt.close()
Quasinormal Fits¶
Each column represents the $\psi_4$ mode as a function of time and rows look at log(amplitude) and phase. The dashed lines correspond to the simulation data where colors indicate the simulation. Solid lines are fit to the data (not quasinormal ringdown fits!).
To Check - If (3,2) mode in BBH is also equally noisy.
In [19]:
#Computing the quasinormal frequencies with given mass and spins using Ringdown fits
import pandas as pd
#Reference - PHYSICAL REVIEW D73,064030 (200
def fit_omega(a_bh, M_bh, f1, f2, f3):
omega_lmn = (f1 + f2*(1.-a_bh)**f3)/M_bh
return omega_lmn
def fit_tau(a_bh, omega_lmn, q1, q2, q3):
tau_lmn = (2./omega_lmn)*(q1+q2*(1.-a_bh)**q3)
return tau_lmn
f1_dict = {'220':1.5251, '210':0.6,'200':0.4437, '2-10':0.3441, '2-20':0.2938,\
'330':1.8956, '320':1.1481,'310':0.8345,'300':0.6873,'3-10':0.5751,'3-20':0.5158,'3-30':0.4673}
f2_dict = {'220':-1.1568, '210':-0.2339,'200':-0.0739, '2-10':0.0293, '2-20':0.0782,\
'330':-1.3043, '320':-0.5552,'310':-0.24205,'300':-0.09282,'3-10':0.02508,'3-20':0.8195,'3-30':0.1296}
f3_dict = {'220':0.1292, '210':0.4175,'200':0.335, '2-10':2.001, '2-20':1.3546,\
'330':0.1818, '320':0.3002,'310':0.4095,'300':0.3479,'3-10':3.136,'3-20':1.4084,'3-30':1.3255}
q1_dict = {'220':0.7, '210':-0.3,'200':4.0, '2-10':2.0, '2-20':1.67,
'330':0.9, '320':0.8313,'310':23.845,'300':6.7841,'3-10':3.0464,'3-20':2.9,'3-30':2.55}
q2_dict = {'220':1.4187, '210':2.3561,'200':-1.955, '2-10':0.1078, '2-20':0.4192,
'330':2.343, '320':2.3773,'310':-20.724,'300':-3.6112,'3-10':0.1162,'3-20':0.3356,'3-30':0.6576}
q3_dict = {'220':-0.499, '210':-0.2277,'200':0.142, '2-10':5.0069, '2-20':1.47,
'330':-0.481, '320':-0.3655,'310':0.03837,'300':0.0948,'3-10':-0.2812,'3-20':2.305,'3-30':1.3378}
#From results obtained below
q2_bbh = ['BBH', 2.0,0.6137, 0.9591]
q2_nsbh = ['NSBH',2.0,0.6846, 0.9672]
q3_bbh = ['BBH', 3.0,0.5405, 0.9712]
q3_nsbh = ['NSBH',3.0,0.564, 0.9746]
q5_bbh = ['BBH', 5.0,0.4166, 0.9824]
q5_nsbh = ['NSBH',5.0,0.4204, 0.9834]
remnant_data = pd.DataFrame(np.array((q2_bbh, q2_nsbh, q3_bbh, q3_nsbh, q5_bbh, q5_nsbh)),columns=('System','q','a','Mh') )
for key in ['210', '220','330']:#,'440']:
f1 = f1_dict[key]
f2 = f2_dict[key]
f3 = f3_dict[key]
q1 = q1_dict[key]
q2 = q2_dict[key]
q3 = q3_dict[key]
a_bh = pd.to_numeric(remnant_data['a'])
M_bh = pd.to_numeric(remnant_data['Mh'])
omega = fit_omega(a_bh, M_bh, f1, f2, f3)
tau = fit_tau(a_bh, omega, q1, q2, q3)
remnant_data['Omega_%s'%key] = omega
remnant_data['tau_%s'%key] = tau
#remnant_data['MOmega_%s'%key] = M_bh*omega
#remnant_data['Tau_by_M%s'%key] = M_bh*1./tau
In [20]:
#QNM Fitting - Thesis Figure
#Psi4 plots: BBH vs NSBH
from collections import defaultdict
i=0
color = ['k','b', 'orangered','m','forestgreen','r', 'b', 'k','coral']
ls = ['-', '--', ':','-.','--',':']
lm_modes = {'2,2':(2,2), '2,1':(2,1),'3,3':(3,3),'3,2':(3,2), '4,4':(4,4)}
y_amp_max = [1e-4, 5e-4, 1e-4, 1e-4]
y_phase_max = [150, 200,250,250]
y_phase_min = [-25, -50,-20,50]
qnm_tau_dict = defaultdict(dict)
qnm_freq_dict = defaultdict(dict)
fig = plt.figure(figsize=(14,20))
nsbh_iter = 0
for lm_idx, lm_key in enumerate(sorted(lm_modes.keys())):
lm=lm_modes[lm_key]
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
col=compare_nsbh[key]['color']#color[i]
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
nsbh_iter+=1
#remnant_data_nsbh = remnant_data[remnant_data['System']=='NSBH']
q = 5
#q_arr = pd.to_numeric(remnant_data_nsbh['q'])
#remnant_data_nsbh = remnant_data_nsbh[q_arr==q]
#omega = np.float(remnant_data_nsbh['Omega_%d%d0'%(lm[0],lm[1])])
#tau = np.float(remnant_data_nsbh['tau_%d%d0'%(lm[0],lm[1])])
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
phase = -1*np.unwrap(np.arctan2(psi4_im22,psi4_re22))
#Compute QNM Modes
#For QNM - A = constant*exp(-t\Ï„) and phase = omega*t+phi_0
#Find the fitting time interval
max_amp_idx = np.argmax(amp==np.amax(amp))
qnm_begin_idx = np.amin(np.where(t22 > t22[max_amp_idx]+30))
qnm_end_idx = np.amin(np.where(t22 > t22[max_amp_idx]+50))
#Define the fitting region
qnm_time_interval = t22[qnm_begin_idx:qnm_end_idx]
qnm_amp_interval = amp[qnm_begin_idx:qnm_end_idx]
qnm_phase_interval = phase[qnm_begin_idx:qnm_end_idx]
#Create the fit for amplitude and phase
logamp_fit_coeff = np.polyfit(qnm_time_interval, np.log(qnm_amp_interval),1)
logamp_fit_data = np.poly1d(logamp_fit_coeff)
#logamp_fit_berti = lambda x: -(1./tau)*x + (np.log(amp)[qnm_begin_idx] + (1./tau)*t22[qnm_begin_idx])
phase_fit_coeff = np.polyfit(qnm_time_interval, qnm_phase_interval,1)
phase_fit_data = np.poly1d(np.array([phase_fit_coeff[0], phase_fit_coeff[1]]))#np.poly1d(phase_fit_coeff)
#phase_fit_berti = lambda x: omega*x + (phase[qnm_begin_idx] - omega*t22[qnm_begin_idx])
#Compute the Quasinormal Frequency and Decay Time constant from the fits
logamp_fit_slope = logamp_fit_coeff[0]
qnm_tau = -1./logamp_fit_slope
qnm_freq = phase_fit_coeff[0]/2./np.pi
qnm_tau_dict[lm_key][labelname] = qnm_tau
qnm_freq_dict[lm_key][labelname] = qnm_freq
#Plot the data and the fits
plt.subplot(5,2,2*lm_idx+1)
amp_plot = plt.semilogy(t22-t22[max_amp_idx], amp, ls='--', c=col, label=labelname)
plt.semilogy(qnm_time_interval-t22[max_amp_idx], np.exp(logamp_fit_data(qnm_time_interval)),ls='-', c=col)
#plt.semilogy(qnm_t_int-t22[max_amp_idx], np.exp(logamp_fit_berti(qnm_t_int)),ls='-', c=color[lm_idx])
if nsbh_iter==3:plt.xlabel(r'$T/M$')
plt.ylabel(r'$|\psi^{ %d,%d}|$'%(lm[0],lm[1]))
plt.xlim(0,130)
plt.ylim(1e-9,4e-4)
plt.legend()
plt.grid(True)
plt.subplot(5,2,2*lm_idx+2)
plt.plot(t22 - t22[max_amp_idx], phase,ls='--', c=col, label=labelname+'(%d,%d)'%(lm[0],lm[1]))
plt.plot(qnm_time_interval - t22[max_amp_idx], phase_fit_data(qnm_time_interval), ls='-', c=col)#, label=labelname+'Fit')
#plt.plot(qnm_t_int - t22[max_amp_idx], phase_fit_berti(qnm_t_int), ls='-', c=color[lm_idx])#, label=labelname+'Fit')
if lm_idx==2:plt.xlabel(r'$T/M$')
plt.ylabel(r'$\phi^{ %d,%d}$'%(lm[0],lm[1]))
plt.xlim(0,130)
#plt.ylim(y_phase_min[nsbh_iter-1], y_phase_max[nsbh_iter-1])
plt.legend(loc=2)
plt.grid(True)
#fig.suptitle('(%s) Mode'%lm_key, fontsize=16, y=1.08)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'QNMFitting_NSBH.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
# QNM Tau & QNM Frequency Comparison - Thesis Figure
color = ['C0','C1','C2','C3','C4','C5']
fig=plt.figure(figsize=(15,6))
plt.subplot(121)
for key1 in sorted(qnm_tau_dict.keys()):
print("\n\nMode: l=%s, m=%s \n"%(key1[0], key1[2]))
print("Simulation \t QNM Tau \t QNM Freqeuncy")
print("-"*50)
for i,key2 in enumerate(sorted(qnm_tau_dict[key1].keys())):
print("%10s \t %g \t %g "%(key2, qnm_tau_dict[key1][key2], qnm_freq_dict[key1][key2]*2*np.pi))
if key1=='2,2':label=key2
else: label= "_nolegend_"
if key2=='Q2-NSBH':continue
plot = plt.scatter(key1, qnm_tau_dict[key1][key2], c=color[i],label=label)
plt.xlabel('(l,m) Mode')
plt.ylabel(r'$\tau_{QNM}$')
plt.grid(True)
plt.legend()
color = ['C0','C1','C2','C3','C4','C5']
plt.subplot(122)
for key1 in sorted(qnm_freq_dict.keys()):
for i,key2 in enumerate(sorted(qnm_freq_dict[key1].keys())):
if key2=='Q2-NSBH':continue
if key1=='2,2':label=key2
else: label= "_nolegend_"
plot = plt.scatter(key1, qnm_freq_dict[key1][key2], c=color[i],label=label)
print('\n\n')
plt.xlabel('(l,m) Mode')
plt.ylabel(r'$f_{QNM}$')
plt.grid(True)
plt.legend()
plt.tight_layout()
#fig.suptitle('Figure 7: QNM Quatities', y=1.08, fontsize=16)
plt.savefig(os.path.join(figdir,'QNMModes.png'), dpi=100)
#plt.show()
plt.close()
Radiated Energy, Angular Momentum and BH Kicks¶
Radiated Quantities¶
Three key observations -
- Radiated energy and angular momentum after merger increase with BH spin.
- At the merger, BH kicks are highest for non-spinning case and gradually decrease for each orientation of the spin. In contrast, post merger kicks are highest for anti-aligned spin.
- BH suffers a strong retardation for very high spin in +z direction which is reflected in decrease in the speed of the black hole. In contrast the speeds increase for the anti-aligned case!
In [21]:
#Radiated Quantities
from scipy.constants import c, G
print('-'*100)
print('Simulation \t Time after merger \t Radiated Energy \t Luminosity \t \tRadiated Jz ')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
t_merge = compare_nsbh[key]['t_merger_src']
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
#filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
if not os.path.exists(filepath):continue
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
time_centered = t - (t_merge + 70)
try:
idx = np.amin(np.where(time_centered>200))
except ValueError as err:
idx=-1
print ('{:<16s} \t {:>4.3g} \t \t {:>4.3g} \t \t{:>4.3g} \t {:>12g} '.format(labelname, time_centered[idx], Erad[idx], Edot[idx], Jz[idx]))
print('\n\n')
print('-'*100)
print('Simulation \t Merger Time \tBH Mass \t Kick (km/s) \t Merger Time+200 BH Mass \t Kick (km/s)')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
t_merge = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
rad_arr = np.array((30,40,50,60,70,75,80,90,100,115,130,145,160,175,190,200))
t_merge = compare_nsbh[key]['t_merger_src']
#filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):continue
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
Kmag = (Kx**2 + Ky**2 + Kz**2)**0.5
time_centered = t - (t_merge+70)
idx_merger = np.amin(np.where(time_centered>0))
idx_postmerger = np.amin(np.where(time_centered>200))
bh_diag_file = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
ihspin_file = os.path.join(datadir, 'ihspin_hn_0.asc')
t_ah, irr_mass = np.loadtxt(bh_diag_file, unpack=True, usecols=(1,26))
t_ih, sx, sy, sz = np.loadtxt(ihspin_file, unpack=True, usecols=(0,1,2,3))
ah_idx_merger = np.amin(np.where(t_ah>t_merge + 0))
ah_idx_postmerger = np.amin(np.where(t_ah>t_merge + 200))
ih_idx_merger = np.amin(np.where(t_ih>t_merge + 0))
ih_idx_postmerger = np.amin(np.where(t_ih>t_merge + 200))
spin_mag_merger = np.linalg.norm(np.array((sx, sy, sz)).T, axis=1)[ih_idx_merger]
hzn_mass_merger = np.sqrt(irr_mass[ah_idx_merger]**2 + 0.25*spin_mag_merger**2/irr_mass[ah_idx_merger]**2)
spin_mag_postmerger = np.linalg.norm(np.array((sx, sy, sz)).T, axis=1)[ih_idx_postmerger]
hzn_mass_postmerger = np.sqrt(irr_mass[ah_idx_postmerger]**2 + 0.25*spin_mag_postmerger**2/irr_mass[ah_idx_postmerger]**2)
#print ('{:<16s} {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.4g} \t {:>1.4g}'.format(labelname,t[idx],Kx[idx]/hzn_mass, Ky[idx]/hzn_mass, \
# Kz[idx]/hzn_mass, hzn_mass, (Kmag[idx]/hzn_mass)*c/1000 ))
print ('{:<16s} {:>1.1f} \t \t {:>1.4g} \t {:>1.3g} \t\t {:>1.1f} \t {:>1.4g} \t {:>1.4g} '.format(labelname,time_centered[idx_merger], \
hzn_mass_merger, \
(Kmag[idx_merger]/hzn_mass_merger)*c/1000,\
time_centered[idx_postmerger],hzn_mass_postmerger, \
(Kmag[idx_postmerger]/hzn_mass_postmerger)*c/1000 ))
Radiated Energy, Luminosity and Angular Momentum¶
System with highest initial BH spin emits maximum energy in gravitational radiation and follows a decreasing trend with spin. In contrast, the luminosity peaks at $a_{BH}=0.25$ and decreases for all other values. This trend differs slightly from $\psi_4^{22}$.
Angular momentum shows two interesting trends - first for non-negative spinning systems, there exist additional angular momentum trapped in the remaining matter. Second, the gravitational radiation is emitted for longer duration in situations with minimal stellar disruptions. In contrast, the loss of angular momentum in gravitational waves stops close to the merger for high spinning cases.
In [22]:
#Energy and Angular Momentum Plot:
i=0
plt.figure(figsize=(15,12))
color = [ 'k','b','orangered','m','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
col = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls='-'
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls='--'
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):
print('File note found for {}'.format(labelname))
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(221)
plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$E_{rad}/M$')
plt.xlim(-400,200)
plt.legend()
plt.subplot(222)
plt.plot(t, Edot,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$\dot{E}_{rad}$')
#plt.xlim(-150,100)
plt.legend()
plt.subplot(223)
plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$J_{rad}/M$')
plt.xlim(-400,200)
plt.legend()
plt.subplot(224)
plt.plot(t, Jzdot,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$\dot{J}_{rad}$')
#plt.xlim(-150,100)
plt.legend()
plt.tight_layout()
#plt.savefig(os.path.join(figdir,'ERad.png'), dpi=100)
plt.show()
plt.close()
In [23]:
#Angular Momenta and BH Spin
import matplotlib.ticker as plticker
i=0
fig, ax = plt.subplots(2,3, figsize=(16,10))
ax = np.concatenate(ax)
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['--', '-', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
mass_ratio = 5
spin_bh = float(labelname.split('=')[1].split('$')[0])
ls = compare_nsbh[key]['ls']
color=compare_nsbh[key]['color']
if spin_bh==-0.5:
ax1 = ax[0]
tmax = 300
#ymin,ymax = -1.5e-3, 1.5e-3
elif spin_bh==0.:
ax1 = ax[1]
tmax = 700
#ymin,ymax = -1.e-3, 1.e-3
elif spin_bh==0.25:
ax1 = ax[2]
tmax = 900
elif spin_bh==0.5:
ax1 = ax[3]
tmax = 1000
elif spin_bh==0.75:
ax1 = ax[4]
tmax = 1200
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?0?.maximum.asc'))[0]
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
tmerge_bbh = compare_nsbh[key]['t_merger_src']
ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?2?.maximum.asc'))[0]
color='b'
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
r=70
#Read Radiated quantities data
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, \
Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
Jz_adm = compare_nsbh[key]['J_adm'][-1]
#M_adm = compare_nsbh[key]['M_adm']
#Read BH Spin data
#t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
R_ex_gwrad=float((os.path.basename(psi4rad_filepath).split('.asc')[0]).split('_r')[1])
t_ih, sz = np.loadtxt(ihspin, unpack=True, usecols=(1,2))
#print("{}: ADM Angular Momentum - {:>2.3g}, Radiated Angular Momentum = {:>2.5g}, BH Spin = {:>2.5g}, \
#BH spin + GW ang mom = {:>2.3g}".format(labelname, Jz_adm, Jz_rad[-1], sz[-1], np.absolute(Jz_rad[-1]) + np.absolute(sz[-1])))
i=0
#print(labelname, np.absolute(J_adm)-np.absolute(Jz_rad[-1]) -np.absolute(sz[-1]))
ax1.plot(t_rad, -np.absolute(Jz_rad )+ np.absolute(Jz_adm),c=color,lw=2,ls='-', label=labelname)
ax1.plot(t_ih, sz, c='k',ls='--', label='BH Spin')
ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$|J^{ADM} - J^{GW}|\,/M^2$', fontsize=18)
ax1.set_xlim(0,tmax)
ax1.legend()
ax1.grid(True)
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
#plt.savefig(os.path.join(figdir, 'Radiated_AngMom.png'), dpi=200)
plt.show()
plt.close()
i=0
ADM Energy and Angular Momentum Conservation¶
We next investigate the conservation of ADM energy and angular momentum to determine the energy and angular momentum trapped in the disk. Few interesting observations -
- While the fraction of ADM energy emitted in gravitational waves increase with BH spin, reverse trend is observed in Angular momentum.
- Due to stellar disruption, the mass and angular momentum of $a_BH>0.5$ continues to grow. In contrast, the remaining systems show a more BBH like behavior.
- The disk reaches upto 2% of ADM energy and 7% of angular momentum for highest spinning case.
In [24]:
#Energy And Angular MomentumPlot:
i=0
plt.figure(figsize=(15,10))
color = [ 'k','b','orangered','m','forestgreen']
ls = ['-', '--', ':','-.','--',':']
print('-'*100)
print('System \t\t M_bh/M_adm \t E_gw/M_adm) \t E_Disk/M_adm \t S_bh/J_adm \t J_gw/J_adm \t J_disk/J_adm')
print('-'*100)
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
col = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls='-'
ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?0?.maximum.asc'))[0]
bhdiag = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ihspin = glob.glob(os.path.join(datadir, 'qlm_mp_j1?2?.maximum.asc'))[0]
bhdiag = os.path.join(datadir, 'BH_diagnostics.ah3.gp')
ls='--'
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):
print('File note found for {}'.format(labelname))
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
J_adm = compare_nsbh[key]['J_adm'][-1]
M_adm = compare_nsbh[key]['M_adm']
bhdiag_filepath = os.path.join(datadir, 'BH_diagnostics.ah1.gp')
t_ah, irr_mass = np.loadtxt(bhdiag, unpack=True, usecols=(1,26))
t_ih, sz = np.loadtxt(ihspin, unpack=True, usecols=(1,2))
spin_splinefit = splrep(t_ih[::8], sz[::8],s=0)
sz_fit = splev(t_ah, spin_splinefit, der=0)
hzn_mass = np.sqrt(irr_mass**2 + 0.25*sz_fit**2/irr_mass**2)
print("{:15s} {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g} \t {:>2.4g}".format(labelname,\
hzn_mass[-1]/M_adm, Erad[-1]/M_adm, 1.0-(hzn_mass[-1]+ Erad[-1])/M_adm, \
sz[-1]/J_adm, Jz_rad[-1]/J_adm, 1.0 - \
(np.absolute(sz[-1]/J_adm) + np.absolute(Jz_rad[-1]/J_adm)) ))
plt.subplot(221)
plt.plot(t-t22[max_amp_idx], Erad/M_adm,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'$E_{rad}/M_{ADM}$')
plt.xlim(-400,200)
plt.legend()
plt.subplot(222)
plt.plot(t_ah-t22[max_amp_idx]+70, hzn_mass/M_adm, ls=ls, c=col, label=labelname )
plt.xlabel(r'$T/M$')
plt.ylabel(r'$M_{BH}/M_{ADM}$')
plt.xlim(-400,200)
plt.legend()
plt.subplot(223)
plt.plot(t-t22[max_amp_idx], Jz_rad/J_adm,ls=ls, c=col, label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'${J}_{rad}/J_{ADM}$')
plt.xlim(-150,100)
plt.legend()
plt.subplot(224)
plt.plot(t_ih-t22[max_amp_idx], sz/J_adm,ls=ls, c=color[i], label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel(r'${S}_{z}/J_{ADM}$')
plt.xlim(-150,100)
plt.ylim(-0.5, 1)
plt.legend()
print('\n')
plt.tight_layout()
#plt.savefig(os.path.join(figdir,'C3_Radiated_Energy_AngMom.png'), dpi=100)
plt.show()
plt.close()
Black Hole Kicks¶
Kicks vary as a function of spin but trends are not very clear. For $a_{BH}=-0.5$, black hole is strongly accelerated and is able to retain its velocity instead of dissipating in gravitational waves. This is also consistent with least energy emitted in gravitational radiation.
As the spin increases, the acceleration is quickly followed by
In [25]:
# BH Kicks
import matplotlib.gridspec as gridspec
import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
fig = plt.figure(constrained_layout=True, figsize=(14,15))
gs = gridspec.GridSpec(ncols=2, nrows=3, figure=fig)
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,1])
ax3 = fig.add_subplot(gs[1,0])
ax4 = fig.add_subplot(gs[1,1])
ax5 = fig.add_subplot(gs[2,0])
ax6 = fig.add_subplot(gs[2,1])
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = '5'
spin = float(labelname.split('=')[1][:-1])
color = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
t_ih3, sx3, sy3, sz3, px3, py3, pz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3,4,5,6))
if spin<0.25:
irrmass_splinefit = splrep(t_ah3[::8],irrmass_ah[::8],s=0)
spinz_splinefit = splrep(t_ih3[::2], sz[::2], s=0)
else :
irrmass_splinefit = splrep(t_ah3[::128],irrmass_ah[::128],s=0)
spinz_splinefit = splrep(t_ih3[::16], sz[::16], s=0)
sz_fit = splev(t_ah3, spinz_splinefit, der=0)
hzn_mass = np.sqrt(irrmass_ah3**2 + 0.25*sz_fit**2/irrmass_ah3**2)
hzn_mass_fit = splrep(t_ah3, hzn_mass, s=0)
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
hzn_mass_evalfit = splev(trad, hzn_mass_fit, der=0)
ax1.plot(trad - (t_merge+70), Kxrad, c=color[i], ls='--')
ax1.plot(t_ih3, px, c=color[i], ls=':')
ax2.semilogy(trad - t_merge-70, Pxdotrad, c=color[i], ls='--')
ax3.plot(trad - (t_merge+70), Kyrad, c=color[i], ls='--')
ax3.plot(t_ih3, py, c=color[i], ls=':')
ax4.semilogy(trad - t_merge-70, Pydotrad, c=color[i], ls='--')
ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2), c=color[i], ls='--')
ax5.plot(t_ih3, np.sqrt(px**2 + py**2 + pz**2), c=color[i], ls=':')
ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color[i], ls='--')
print(labelname,bbh_irrmass_bh3, bbh_spinmag_bh3[-1])
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx, sy, sz, px, py, pz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3,4,5,6))
if spin<0.25:
irrmass_splinefit = splrep(t_ah[::8],irrmass_ah[::8],s=0)
spinz_splinefit = splrep(t_ih[::2], sz[::2], s=0)
else :
irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
spinz_splinefit = splrep(t_ih[::16], sz[::16], s=0)
sz_fit = splev(t_ah, spinz_splinefit, der=0)
hzn_mass = np.sqrt(irrmass_ah**2 + 0.25*sz_fit**2/irrmass_ah**2)
hzn_mass_fit = splrep(t_ah, hzn_mass, s=0)
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
hzn_mass_evalfit = splev(trad-70, hzn_mass_fit, der=0)
ax1.plot(trad - (t_merge+70), Kxrad/hzn_mass_evalfit, c=color, ls='-',label=labelname)
#ax1.plot(t_ih - t_merge, px, c=color, ls=':')
ax2.plot(trad - t_merge-70, Pxdotrad, c=color, ls='-',label=labelname)
ax3.plot(trad - (t_merge+70), Kyrad/hzn_mass_evalfit, c=color, ls='-',label=labelname)
#ax3.plot(t_ih - t_merge, py, c=color, ls=':')
ax4.plot(trad - t_merge-70, Pydotrad, c=color, ls='-',label=labelname)
ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)/hzn_mass_evalfit, c=color, ls='-',label=labelname)
#ax5.plot(t_ih - t_merge, np.sqrt(px**2 + py**2 + pz**2), c=color, ls=':')
ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color, ls='-',label=labelname)
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ K_x$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{K}_x$', fontsize=20)
ax3.set_xlabel(r'$T/M$', fontsize=20)
ax3.set_ylabel(r'$K_y$', fontsize=20)
ax4.set_xlabel(r'$T/M$', fontsize=20)
ax4.set_ylabel(r'$\dot{K}_y$', fontsize=20)
ax5.set_xlabel(r'$T/M$', fontsize=20)
ax5.set_ylabel(r'$K/M$', fontsize=20)
ax6.set_xlabel(r'$T/M$', fontsize=20)
ax6.set_ylabel(r'$\dot{K}$', fontsize=20)
ax1.set_xlim(-50,80)
ax2.set_xlim(-50,80)
ax3.set_xlim(-50,80)
ax4.set_xlim(-50,80)
ax5.set_xlim(-50,80)
ax6.set_xlim(-50,80)
#ax2.set_ylim(2e-8,2e-4)
#ax4.set_ylim(2e-8,2e-4)
ax6.set_ylim(2e-8,2e-4)
ax1.legend(fontsize=12, loc=4)
ax2.legend(fontsize=12, loc=1)
ax3.legend(fontsize=12, loc=4)
ax4.legend(fontsize=12)
ax5.legend(fontsize=12)
ax6.legend(fontsize=12)
ax1.grid(True)
ax2.grid(True)
ax3.grid(True)
ax4.grid(True)
ax5.grid(True)
ax6.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_Kick.png'), dpi=100)
plt.show()
plt.close()
QuasiLocal Quantities¶
Black Hole and Neutron Star Mass¶
In [26]:
#BH Horizon Mass Plots - Comparison of mass and spin of final black hole at 300M after merger
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_irrmass_dict = {}
bbh_hznmass_dict = {}
print("System \t \t Time \t Init. Irr Mass Final Irr Mass \t Final Spin(a) \tFinal Hzn Mass \n ")
print('-'*115)
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = '5'
spin = float(labelname.split('=')[1][:-1])
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
bbh_irrmass_bh1 = irrmass_ah1[0]
bbh_spinmag_bh1 = (sx1**2 + sy1**2 + sz1**2)**0.5
bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
#Set all the simulations to same final time before comparing mass and spins
area_ah3 = area_ah3 [t_ah3 < t_comp]
irrmass_ah3 = irrmass_ah3 [t_ah3 < t_comp]
arealrad_ah3 = arealrad_ah3[t_ah3 < t_comp]
t_ah3 = t_ah3 [t_ah3 < t_comp]
sx3 = sx3 [t_ih3 < t_comp]
sy3 = sy3 [t_ih3 < t_comp]
sz3 = sz3 [t_ih3 < t_comp]
t_ih3 = t_ih3[t_ih3 < t_comp]
bbh_irrmass_bh3 = irrmass_ah3[-1]
bbh_spinmag_bh3 = (sx3**2 + sy3**2 + sz3**2)**0.5
bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
bbh_irrmass_dict[massratio] = bbh_irrmass_bh3
bbh_hznmass_dict[massratio] = bbh_hznmass_bh3
print("%-10s \t %g \t %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
bbh_hznmass_bh3, bbh_hznmass_bh3 - bbh_hznmass_bh3))
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 400.
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
mass_growth_rate = np.gradient(irrmass_ah[::8])/np.gradient(t_ah[::8])
#Set all the simulations to same final time before comparing mass and spins
area_ah = area_ah [t_ah<t_comp]
irrmass_ah = irrmass_ah [t_ah<t_comp]
arealrad_ah = arealrad_ah[t_ah<t_comp]
t_ah = t_ah [t_ah<t_comp]
sx = sx [t_ih<t_comp]
sy = sy [t_ih<t_comp]
sz = sz [t_ih<t_comp]
t_ih = t_ih[t_ih<t_comp]
nsbh_irrmass_bh1 = irrmass_ah[0]
nsbh_spinmag_bh1 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
nsbh_irrmass_bh3 = irrmass_ah[-1]
nsbh_spinmag_bh3 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
print("%-1s \t %.1f \t %g \t %g \t \t%g \t %g "%(labelname[1:-1], t_ah[-1]-t_merge, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
nsbh_dimspin_bh3, nsbh_hznmass_bh3))
In [42]:
#Masss andd Spin Plot (Paper)
import matplotlib.gridspec as gridspec
import matplotlib as mpl
mpl.rc('lines', linewidth=2, color='r')
mpl.rc('font', size=18)
mpl.rc('axes', labelsize=18, grid=True)
mpl.rc('xtick', labelsize=14)
mpl.rc('ytick', labelsize=14)
mpl.rc('legend', fontsize=16)
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
fig = plt.figure(constrained_layout=True, figsize=(14,10))
gs = gridspec.GridSpec(ncols=2, nrows=2, figure=fig)
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,1])
ax3 = fig.add_subplot(gs[1,0])
ax4 = fig.add_subplot(gs[1,1])
#ax5 = fig.add_subplot(gs[2,0])
#ax6 = fig.add_subplot(gs[2,1])
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = '5'
spin = float(labelname.split('=')[1][:-1])
color = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
t_ih3, sx3, sy3, sz3,px3,py3,pz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3,4,5,6))
#Set all the simulations to same final time before comparing mass and spins
irrmass_ah3 = irrmass_ah3 [t_ah3 < t_comp]
t_ah3 = t_ah3 [t_ah3 < t_comp]
sx3 = sx3 [t_ih3 < t_comp]
sy3 = sy3 [t_ih3 < t_comp]
sz3 = sz3 [t_ih3 < t_comp]
t_ih3 = t_ih3[t_ih3 < t_comp]
bbh_irrmass_bh3 = irrmass_ah3[-1]
bbh_spinmag_bh3 = (sx3**2 + sy3**2 + sz3**2)**0.5
bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
bbh_dimless_spin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
ax1.axhline(y=bbh_irrmass_bh3,c=color[i], ls='--')
ax3.axhline(y=bbh_spinmag_bh3[-1],c=color[i], ls='--')
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
#ax5.plot(trad - (t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2), c=color[i], ls='--')
#ax5.plot(t_ih3, np.sqrt(px**2 + py**2 + pz**2), c=color[i], ls=':')
#ax6.semilogy(trad - t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color[i], ls='--')
print(labelname,bbh_irrmass_bh3, bbh_spinmag_bh3[-1])
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
if spin<0.25:
irrmass_splinefit = splrep(t_ah[::8],irrmass_ah[::8],s=0)
spinz_splinefit = splrep(t_ih[::2], sz[::2], s=0)
else :
irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
spinz_splinefit = splrep(t_ih[::16], sz[::16], s=0)
#Compute mass growth rate on every 8th iteration to reduce the high frequency noise in the derivatives
#Truncation and round off errors leads to high frequency noise which gets amplified on computing the derivatives
#Computing the derivative on a larger time step averages out these effect to zero)
mass_growth_rate_fit = splev(t_ah, irrmass_splinefit, der=1)
mass_growth_rate = np.gradient(irrmass_ah[::32])/np.gradient(t_ah[::32])
spin_growth_rate_fit = splev(t_ih, spinz_splinefit, der=1)
spin_growth_rate = np.gradient(sz)/np.gradient(t_ih)
sz_fit = splev(t_ah, spinz_splinefit, der=0)
hzn_mass = np.sqrt(irrmass_ah**2 + 0.25*sz_fit**2/irrmass_ah**2)
hzn_mass_fit = splrep(t_ah, hzn_mass, s=0)
ax1.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color, ls='-',label=labelname)
#ax1.plot(t_ah-t_merge, splev(t_ah, irrmass_splinefit), color[i], ls='--')
ax2.semilogy((t_ah-t_merge), mass_growth_rate_fit, c=color, ls='-',label=labelname)
#ax2.semilogy((t_ah-t_merge)[::32], mass_growth_rate, c=color[i], ls='--',label=labelname[:2])
ax3.plot(t_ih-t_merge, sz-sz[0], c=color, ls='-',label=labelname)
ax4.semilogy(t_ih-t_merge, spin_growth_rate, c=color, ls='-',lw=2,label=labelname)
#ax4.semilogy(t_ih-t_merge, spin_growth_rate_fit, c=color[i-3], ls='--', lw=2,label=labelname[:2])
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
hzn_mass_evalfit = splev(trad, hzn_mass_fit, der=0)
#ax5.plot(trad-(t_merge+70), np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)/hzn_mass_evalfit, c=color, ls='-',label=labelname)
#ax5.semilogy(trad, np.sqrt(Px**2 + Py**2 + Pz**2), c=color[i], ls='--',label=labelname[:2])
#ax6.semilogy(trad-t_merge-70, np.sqrt(Pxdotrad**2 + Pydotrad**2 + Pzdotrad**2), c=color, ls='-',label=labelname)
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ \Delta M_i/M$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{M}_i$', fontsize=20)
ax3.set_xlabel(r'$T/M$', fontsize=20)
ax3.set_ylabel(r'$\Delta S_z/M^2$', fontsize=20)
ax4.set_xlabel(r'$T/M$', fontsize=20)
ax4.set_ylabel(r'$\dot{S}_z/M$', fontsize=20)
#ax5.set_xlabel(r'$T/M$', fontsize=20)
#ax5.set_ylabel(r'$|P/M|$', fontsize=20)
#ax6.set_xlabel(r'$T/M$', fontsize=20)
#ax6.set_ylabel(r'$|\dot{P}|$', fontsize=20)
ax1.set_xlim(-200,300)
ax2.set_xlim(-50,80)
ax3.set_xlim(-200,300)
ax4.set_xlim(-50,80)
#ax5.set_xlim(-200,300)
#ax6.set_xlim(-50,80)
#ax1.set_ylim(-0.01,0.005)
ax2.set_ylim(2e-5,0.1)
ax4.set_ylim(2e-5,0.2)
#ax6.set_ylim(2e-8,2e-4)
ax1.legend(loc=4)
ax2.legend()
ax3.legend()
ax4.legend()
#ax5.legend()
#ax6.legend()
#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)
#ax5.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax1.grid(True)
ax2.grid(True)
ax3.grid(True)
ax4.grid(True)
#ax5.grid(True)
#ax6.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
plt.show()
plt.close()
Rest Mass outside BH¶
In [31]:
#Rest Mass outside the horizon
#subplt_num=1
print("System \t\t Time - T_merger\t Mass (R>1.5)/M_* ")
print('-'*100)
plt.figure(figsize=(8,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = '5'
spin = float(labelname.split('=')[1][:-1])
t_merge = compare_nsbh[key]['t_merger_src']
color = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='BBH'):
continue
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
#Volume Integrals Data
vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
time, m_ns, m_bh3, m_0d7, m_1d4, m_20, m_35 = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5,6,7,8,9,10))
#Choose merger point as reference for calculation
time_centered = time-t_merge
idx = np.where(time_centered<400)
time_centered = time_centered[idx]
if compare_nsbh[key]['symmetry']:
m_ns = 2.*m_ns[idx] #accounting for symmetries
m_bh3 = 2.*m_bh3[idx] #accounting for symmetries
m_0d7 = 2.*m_0d7[idx] #accounting for symmetries
m_1d4 = 2.*m_1d4[idx] #accounting for symmetries
m_20 = 2.*m_20[idx] #accounting for symmetries
m_35 = 2.*m_35[idx] #accounting for symmetries
#This is approximate mass and we are ignoring the mass in sphere of 0.7M due to lack of data. This mass
#will be of order of 1e-15 and hence can be neglected
restmass_atm = m_0d7[0] - m_ns[0]
#print('Atmosphere rest mass = %g (M) = %g (M_*), m_ns=%g'%(restmass_atm, restmass_atm/m_ns[0], m_ns[0]))
mass_outside_hzn = m_1d4[idx] - restmass_atm
mass_outside_R35 = m_35[idx] - m_35[0]
mass_outside_R20 = m_20[idx] - m_20[0]
m_disk_R35_M = mass_outside_hzn - mass_outside_R35
m_disk_R20_M = mass_outside_hzn - mass_outside_R20
m_disk_R35_Mstar = m_disk_R35_M/m_ns[0]
m_disk_R20_Mstar = m_disk_R20_M/m_ns[0]
m_ejecta_R35_M = mass_outside_R35
mass_outhzn_Mstar = mass_outside_hzn/m_ns[0]
print("%-10s \t %.2f \t \t %.2g "%(labelname[1:-1], time_centered[-1], mass_outhzn_Mstar[-1]))
plt.semilogy(time_centered, mass_outhzn_Mstar, c=color, ls='-', label=labelname)
#plt.semilogy(time_centered, m_disk_R35_Mstar, c=color[i-3], ls='--', label=labelname[:2]+' (R = 35M)')
plt.xlabel(r'$T/M$')
plt.ylabel(r'$M_{td}/M_o$')
plt.ylim(1e-5,2)
plt.xlim(-100,+400)
#plt.tick_params(left=True, right=True, which='both')
#plt.axhline(y=0.06)
plt.legend()
plt.grid(True)
#plt.savefig(os.path.join(figdir, 'ReviewFigures/RevFig_RestMassOutsideBH.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
#Foucart Fits. -
# arXiv:1807.00011 - Equation 4
alpha_fit = 0.406
beta_fit = 0.139
gamma_fit = 0.255
delta_fit = 1.761
def remnant_mass_fit(nsbh_dict):
comp = nsbh_dict['comp']
q = nsbh_dict['q']
a_bh = nsbh_dict['a_bh']
eta = q/(1. + q)**2.
z1 = 1. + ((1. - a_bh**2)**(1./3))*((1. + a_bh)**(1./3) + (1. - a_bh)**(1./3))
z2 = np.sqrt(3.*a_bh**2. + z1**2.)
isco_radius_hat = 3. + z2 - np.sign(a_bh)*np.sqrt((3. - z1)*(3. + z1 + 2*z2))
remnant_mass = (max(alpha_fit*(1-2.0*comp)/(eta**(1./3)) - beta_fit*isco_radius_hat*comp/eta + gamma_fit,\
0))**delta_fit
return remnant_mass
q2_dict = {'comp':0.1529, 'q':2, 'a_bh':0.} #0.683}
q3_dict = {'comp':0.1536, 'q':3, 'a_bh':0.} #0.563}
q5_dict = {'comp':0.1543, 'q':5, 'a_bh':0.} #0.4203}
q7_dict = {'comp':0.145 , 'q':5, 'a_bh':0.0}
print('Q2: Remnant mass = %g'%remnant_mass_fit(q2_dict))
print('Q3: Remnant mass = %g'%remnant_mass_fit(q3_dict))
print('Q5: Remnant mass = %g'%remnant_mass_fit(q5_dict))
print('Q7 (Test): Remnant mass = %g'%remnant_mass_fit(q7_dict))
Multipole Moments¶
Here, we look at the multipole moments. The $M^i_{mp}$ refers to mass multipoles while $J^i_{mp}$ refers to spin multipoles. For mass multipoles, the odd moments are 0 while for spin multipoles, the even moments are zero. The zeroeth mass multipole is basically the Christodoulou mass, as can be seen in the first plot and zeroeth spin multipole in the spin angular momentum of black hole.
Things to look into
- Comparison between BBH and NSBH mass multipoles.
- Does accreting matter affect the multipole moments of black hole?
- Why do higher moments continue to exist (i.e. nonzero) for final black hole?
- Read - Isolated and dynamical horizons and their applications, Section 5.3.
In [46]:
#Quasi Local Moments
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_ah3_mass = 0
#print("System \t\t\t Initial BH Mass \t Final BH Mass \t NSBH - BBH \n ")
fig1, ax1 = plt.subplots(2,2,figsize=(15,8))
ax1 = np.concatenate(ax1)
fig2, ax2 = plt.subplots(2,2,figsize=(15,8))
ax2 = np.concatenate(ax2)
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
spin = float(labelname.split('=')[1][:-1])
t_merge = compare_nsbh[key]['t_merger_src']
color = compare_nsbh[key]['color']
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
qlm_mass_files = sorted(glob.glob(os.path.join(datadir, 'qlm_mp_m*.norm2.asc')))
if len(qlm_mass_files)>0:
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
t_ah = np.around(t_ah, decimals=4)
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx,sy,sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
t_ih = np.around(t_ih, decimals=4)
t_common = np.intersect1d(t_ah, t_ih)
ahf_idx = np.isin(t_ah,t_common)
ih_idx = np.isin(t_ih,t_common)
spinmag = np.sqrt(sx**2 + sy**2 + sz**2)
hznmass_ah = (irrmass_ah[ahf_idx]**2 + 0.25*spinmag[ih_idx]**2/irrmass_ah[ahf_idx]**2)**0.5
for qlm_file in qlm_mass_files:
mass_mp_label = int((os.path.basename(qlm_file).split('[')[0])[-1])
t, mass_moment = np.loadtxt(qlm_file, unpack=True, usecols=(1,2))
if mass_mp_label%2==1 or mass_mp_label>7: continue
ax = ax1[mass_mp_label//2]
ax.plot(t-t_merge,mass_moment, c=color, label='%s'%(labelname))
ax.set_xlabel('Time')
ax.set_ylabel(r'$M_{}$'.format(mass_mp_label))
ax.legend()
ax.grid(True)
#ax.set_ylim(-0.4,0.4)
ax.set_xlim(-50,60)
#if len(qlm_mass_files)>0: ax11.plot(t_ah[ahf_idx], hznmass_ah, ls='--', c='k',label='%s:BH Mass'%(labelname))
qlm_am_files = sorted(glob.glob(os.path.join(datadir, 'qlm_mp_j*.norm2.asc')))
for qlm_file in qlm_am_files:
am_mp_label = int((os.path.basename(qlm_file).split('[')[0])[-1])
if am_mp_label%2==0: continue
t, j_moment = np.loadtxt(qlm_file, unpack=True, usecols=(1,2))
ax = ax2[(am_mp_label-1)//2]
ax.plot(t-t_merge, j_moment, c=color,label='%s'%(labelname))
ax.set_xlabel('Time')
ax.set_ylabel(r'$J_{}$'.format(am_mp_label))
ax.legend()
ax.grid(True)
ax.set_ylim(-0.3,1)
ax.set_xlim(-50,70)
fig1.tight_layout()
fig1.show()
#plt.close()
fig2.tight_layout()
fig2.show()
#plt.close()
Other resources -¶
Additional Plots¶
In [ ]:
simdir1 = '/nethome/numrel/datafiles/Waveforms/HR-series/q5.00/D10_q5.00_a-0.71_0.00_m240/'
psi4_filepath = os.path.join(simdir1, 'mp_WeylScal4::Psi4i_l2_m2_r70.00.asc')
t_psi4, re_psi4, im_psi4 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(re_psi4**2 + im_psi4**2)
t_max_amp = t_psi4[amp == np.amax(amp)]
psi4rad_filepath = os.path.join(simdir1, 'psi4analysis_r70.00.asc') #'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
Kmag = np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)
plt.plot(trad - t_max_amp, Kmag)
simdir2 = '/nethome/numrel/datafiles/Waveforms/HR-series/q5.00/D10_q5.00_a0.6_0.0_m240/'
psi4_filepath = os.path.join(simdir2, 'mp_WeylScal4::Psi4i_l2_m2_r70.00.asc')
t_psi4, re_psi4, im_psi4 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(re_psi4**2 + im_psi4**2)
t_max_amp = t_psi4[amp == np.amax(amp)]
psi4rad_filepath = os.path.join(simdir2, 'psi4analysis_r70.00.asc') #'psi4radAway_l6_r70.00.asc')
trad, Erad, Kxrad, Kyrad, Kzrad, Edotrad, Pxdotrad, Pydotrad, Pzdotrad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8))
Kmag = np.sqrt(Kxrad**2 + Kyrad**2 + Kzrad**2)
plt.plot(trad - t_max_amp, Kmag)
plt.xlim(-50,100)
plt.show()
Orbital Frequency¶
In [ ]:
# Orbital Frequency and rate of change of separation -
plt.figure(figsize=(14,8))
print("\nNote - Mixed binary derivatives are extremely noise. Hence, to tackle the noise, we use spline interpolation \
to fit the original data and compute the derivatives using these fitting functions. While they do provide a \
good approximation to original fits (Tested with BBH), they fail to capture the high frequency oscillations in rdot and omega \
for mixed binaries")
print('-------------'*9)
print('\n\n')
for i, key in enumerate(sorted(compare_nsbh.keys())):
t = orbdyn_dict[key]['time']
label = compare_nsbh[key]['label']
#plt.plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['rdot'], ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
plt.plot(orbdyn_dict[key]['time'], splev(orbdyn_dict[key]['time'], orbdyn_dict[key]['sepfit'], der=1), ls=compare_nsbh[key]['ls'], c=compare_nsbh[key]['color'],label=label+'Spline Fit')
plt.xlabel('Time', fontsize=16)
plt.ylabel('dr/dt', fontsize=16)
plt.legend()
plt.title(r'Figure 3 - $dR_{sep}/dt$', fontsize=20)
plt.xlim(0,900)
plt.ylim(-0.025,0.015)
plt.show()
plt.close()
print('\n\n')
fig = plt.figure(figsize=(15,7))
for i,key in enumerate(sorted(compare_nsbh.keys())):
label = compare_nsbh[key]['label']
omega_fit = splev(orbdyn_dict[key]['time'], orbdyn_dict[key]['phasefit'], der=1)
if 'Q2-NSBH' in label:
time = orbdyn_dict[key]['time'][:-140]
omega_fit = omega_fit[:-140]
else:
time = orbdyn_dict[key]['time']
#print("%s: Initial Orbital Frequency = %g"%(label, orbdyn_dict[key]['omega'][0]))
plt.subplot(121)
#plt.plot(orbdyn_dict[key]['time'],orbdyn_dict[key]['omega'],ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
plt.plot(time, omega_fit, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label+' Spline Derivative')
plt.xlabel('Time/M', fontsize=16)
plt.ylabel(r'$M\Omega_{orb}$', fontsize=16)
plt.ylim(0.022,0.4)
plt.legend()
plt.subplot(122)
#plt.plot(orbdyn_dict[key]['time'], orbdyn_dict[key]['omega']/2./np.pi, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label)
plt.plot(time, omega_fit, ls=compare_nsbh[key]['ls'],c=compare_nsbh[key]['color'],label=label+' Spline Derivative')
plt.xlabel('Time/M', fontsize=16)
plt.ylabel(r'$M\Omega_{orb}$ ', fontsize=16)
plt.legend()
plt.xlim(0,600)
#plt.ylim(0.022,0.1)
plt.tight_layout()
fig.suptitle('Figure 4 - Orbital Frequency', y=1.08, fontsize=20)
#plt.show()
plt.close()
Gravitational Waves¶
In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple
i=0
nsbh_color = ['k','b','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','orangered','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']
fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if 'Q2' not in labelname: continue
if 'Q2' in labelname:
q=2
tmax = 900
l = 0.242
b = 0.892
w = 0.26
h = 0.065
if 'Q3' in labelname:
q=3
tmax = 1000
l = 0.228
b = 0.89
w = 0.235
h = 0.067
if 'Q5' in labelname:
q=5
tmax = 1200
l = 0.203
w = 0.213
h = 0.0706
b = 0.224
ls = compare_nsbh[key]['ls']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
color=nsbh_color
ls='-'
lw=1.5
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
tmerge_bbh = compare_nsbh[key]['t_merger_src']
color=bbh_color
ls='--'
lw=2
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
#Read Psi4 data
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
#Read Radiated quantities data
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
#Read BH Spin data
t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
i=0
ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
#ax1.plot(t22,amp,ls=ls,c=color)
#ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
ax1.set_xlim(0,tmax)
ax1.set_ylim(-6e-2,6e-2)
ax1.legend(loc=3)
ax1.grid(True)
ax1.axvline(150, ls='--', c='grey')
ax1.axvline(451, ls='--', c='grey')
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax1.xaxis.set_major_locator(loc)
i+=1
ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
if 'NSBH' in labelname:
labelplot1 = ax2.plot([],[],c=color[i], ls=ls, label=r'$J_z$')
labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
ax2.set_xlabel(r'$T/M$')
if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
ax2.set_xlim(0,tmax)
ax2.grid(True)
ax2.legend(loc=3)
#ax2.axvline(150, ls='--', c='grey')
#ax2.axvline(450, ls='--', c='grey')
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax2.xaxis.set_major_locator(loc)
#Twin plot in second plot
i+=1
ax21 = ax2.twinx() # instantiate a second axes that shares the same x-axis
if 'NSBH' in labelname:
Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
#ax21.yaxis.set_ticklabels([])
ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i]) # we already handled the x-label with ax1
#ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
else:
print(sz_bh[-1])
Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
#ax21.yaxis.set_ticklabels([])
ax21.tick_params(axis='y', labelright='off',right=False)
ax21.minorticks_off()
ax21.set_ylim(-0.05,0.7)
#ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
#
loc = plticker.MultipleLocator(base=0.15) # this locator puts ticks at regular intervals
ax21.yaxis.set_major_locator(loc)
plt.grid(False)
#Inset in Plot 1
i=0
left, bottom, width, height = [l, b, w, h]
ax12 = fig.add_axes([left, bottom, width, height])
min_idx = np.amin(np.argwhere(t22>=150))
max_idx = np.amin(np.argwhere(t22>=451))
ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
ax12.set_xlim(150, 450)
ax12.set_ylim(-1.3e-3, 1.3e-3)
ax12.yaxis.tick_right()
#ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
ax12.grid(b=True, which='major', axis='x')
ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax12.tick_params(axis='y', labelsize=11)
ax12.yaxis.offsetText.set_fontsize(11)
ax12.minorticks_off()
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.xticks([])
#plt.yticks([])
plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple
i=0
nsbh_color = ['k','b','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','orangered','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']
fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if 'Q3' not in labelname: continue
if 'Q2' in labelname:
q=2
tmax = 900
l = 0.242
b = 0.883
w = 0.26
h = 0.065
if 'Q3' in labelname:
q=3
tmax = 1000
l = 0.228
b = 0.89
w = 0.235
h = 0.067
if 'Q5' in labelname:
q=5
tmax = 1200
l = 0.203
w = 0.213
h = 0.0706
b = 0.224
ls = compare_nsbh[key]['ls']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
color=nsbh_color
ls='-'
lw=1.5
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
tmerge_bbh = compare_nsbh[key]['t_merger_src']
color=bbh_color
ls='--'
lw=2
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
#Read Psi4 data
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
#Read Radiated quantities data
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
#Read BH Spin data
t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
i=0
ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
#ax1.plot(t22,amp,ls=ls,c=color)
#ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
ax1.set_xlim(0,tmax)
ax1.legend(loc=3)
ax1.grid(True)
ax1.axvline(150, ls='--', c='grey')
ax1.axvline(451, ls='--', c='grey')
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax1.xaxis.set_major_locator(loc)
i+=1
ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
if 'NSBH' in labelname:
labelplot1 = ax2.plot([],[],c=color[i], ls=ls, label=r'$J_z$')
labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
ax2.set_xlabel(r'$T/M$')
if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
ax2.set_xlim(0,tmax)
ax2.grid(True)
ax2.legend(loc=3)
#ax2.axvline(150, ls='--', c='grey')
#ax2.axvline(450, ls='--', c='grey')
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax2.xaxis.set_major_locator(loc)
#Twin plot in second plot
i+=1
ax21 = ax2.twinx() # instantiate a second axes that shares the same x-axis
if 'NSBH' in labelname:
Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
#ax21.yaxis.set_ticklabels([])
ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i]) # we already handled the x-label with ax1
#ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
else:
print(sz_bh[-1])
Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
#ax21.yaxis.set_ticklabels([])
ax21.tick_params(axis='y', labelright='off',right=False)
ax21.minorticks_off()
ax21.set_ylim(-0.05,0.6)
#ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
#
loc = plticker.MultipleLocator(base=0.15) # this locator puts ticks at regular intervals
ax21.yaxis.set_major_locator(loc)
plt.grid(False)
#Inset in Plot 1
i=0
left, bottom, width, height = [l, b, w, h]
ax12 = fig.add_axes([left, bottom, width, height])
min_idx = np.amin(np.argwhere(t22>=150))
max_idx = np.amin(np.argwhere(t22>=451))
ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
ax12.set_xlim(150, 450)
ax12.set_ylim(-1e-3, 1e-3)
ax12.yaxis.tick_right()
#ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
ax12.grid(b=True, which='major', axis='x')
ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax12.tick_params(axis='y', labelsize=11)
ax12.yaxis.offsetText.set_fontsize(11)
ax12.minorticks_off()
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.xticks([])
#plt.yticks([])
plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
In [ ]:
#Re(Psi4) and Amplitude plots: Thesis Figures
import matplotlib.ticker as plticker
from matplotlib.legend_handler import HandlerLine2D, HandlerTuple
i=0
nsbh_color = ['k','mediumblue','r','darkviolet','cadetblue', 'r','k','m','darkolivegreen','coral','forestgreen', 'blueviolet']
bbh_color = ['dimgray','royalblue','tomato','darkorchid', 'darkcyan', 'teal','orangered', 'grey',]
ls = ['--', '-', ':','-.','--',':']
fig, (ax1,ax2) = plt.subplots(2, 1, figsize=(10,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if 'Q5' not in labelname: continue
if 'Q2' in labelname:
q=2
tmax = 900
l = 0.242
b = 0.883
w = 0.26
h = 0.065
if 'Q3' in labelname:
q=3
tmax = 1000
l = 0.228
b = 0.89
w = 0.235
h = 0.067
if 'Q5' in labelname:
q=5
tmax = 1200
l = 0.2105
b = 0.8916
w = 0.1983
h = 0.065
ls = compare_nsbh[key]['ls']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
color=nsbh_color
ls='-'
lw=1.5
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
tmerge_bbh = compare_nsbh[key]['t_merger_src']
color=bbh_color
ls='--'
lw=2
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_2.asc')
#Read Psi4 data
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
R_ex_gw=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
#Read Radiated quantities data
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t_rad, Erad_rad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad, Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
#Read BH Spin data
t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
R_ex_gwrad=float((os.path.basename(psi4_filepath).split('.asc')[0]).split('_r')[1])
i=0
ax1.plot(t22, r*psi4_re22,c=color[i],ls=ls, lw=lw)#, label=labelname)
labelplot = ax1.plot([],[],c='k', ls=ls, label=labelname[3:])
#ax1.plot(t22,amp,ls=ls,c=color)
#ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$Re(\psi_4)\,r\,M$')
ax1.set_xlim(0,tmax)
ax1.set_ylim(-3e-2, 3e-2)
ax1.legend(loc=3)
ax1.grid(True)
ax1.axvline(150, ls='--', c='grey')
ax1.axvline(450, ls='--', c='grey')
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax1.xaxis.set_major_locator(loc)
i+=1
ax2.plot(t_rad, Jz_rad,c=color[i],ls=ls, lw=lw)
if 'NSBH' in labelname:
labelplot1 = ax2.plot([],[],c=color[i], ls=ls, label=r'$J_z$')
labelplot2 = ax2.plot([],[],c=color[i+1], ls=ls, label=r'$S_z$')
ax2.set_xlabel(r'$T/M$')
if 'NSBH' in labelname: ax2.set_ylabel(r'$J_z/M^2$')#,color=color[i])
ax2.set_xlim(0,tmax)
ax2.grid(True)
ax2.legend(loc=3)
#ax2.axvline(150, ls='--', c='grey')
#ax2.axvline(450, ls='--', c='grey')
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
loc = plticker.MultipleLocator(base=100) # this locator puts ticks at regular intervals
ax2.xaxis.set_major_locator(loc)
#Twin plot in second plot
i+=1
ax21 = ax2.twinx() # instantiate a second axes that shares the same x-axis
if 'NSBH' in labelname:
Sz_plot = ax21.plot(t_bh, sz_bh, lw=lw, color=color[i], label=r'$S_z$')
#ax21.yaxis.set_ticklabels([])
ax21.set_ylabel(r'$S_z/M^2$')#, color=color[i]) # we already handled the x-label with ax1
#ax21.legend()#, numpoints=1,handler_map={tuple: HandlerTuple(ndivide=None)})
else:
print(sz_bh[-1])
Sz_plot = ax21.axhline(y = sz_bh[-1], c=color[i], lw=lw, ls=ls)
#ax21.yaxis.set_ticklabels([])
ax21.tick_params(axis='y', labelright='off',right=False)
ax21.minorticks_off()
ax21.set_ylim(-0.05,0.5)
#ax21.yaxis.set_ticklabels(np.around(np.arange(0, 0.62,0.15), decimals=2))
#
loc = plticker.MultipleLocator(base=0.1) # this locator puts ticks at regular intervals
ax21.yaxis.set_major_locator(loc)
plt.grid(False)
#Inset in Plot 1
i=0
left, bottom, width, height = [l, b, w, h]
ax12 = fig.add_axes([left, bottom, width, height])
min_idx = np.amin(np.argwhere(t22>=150))
max_idx = np.amin(np.argwhere(t22>=451))
ax12.plot(t22, psi4_re22*r,c=color[i],ls=ls, lw=lw)
ax12.axhline(y=0,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=200,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=300,c='grey', alpha=0.5, lw=1)
ax12.axvline(x=400,c='grey', alpha=0.5, lw=1)
ax12.set_xlim(150, 450)
ax12.set_ylim(-1e-3, 1e-3)
ax12.yaxis.tick_right()
#ax12.tick_params(axis='both',labelbottom='off', labelright='off',bottom=False,right=False)
ax12.grid(b=True, which='major', axis='x')
ax12.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax12.tick_params(axis='y', labelsize=11)
ax12.yaxis.offsetText.set_fontsize(11)
ax12.minorticks_off()
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.xticks([])
#plt.yticks([])
plt.tight_layout()
#fig.suptitle('Figure 1: BBH vs NSBH', fontsize=20, y=1.022)
print(os.path.join(figdir, 'Psi4Plots_q%d.png'%q))
plt.savefig(os.path.join(figdir, 'Psi4Plots_q%d.png'%q), dpi=200)
plt.show()
plt.close()
In [ ]:
# Q2 psi4 frequency oscillations
from scipy import signal
plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
q = float((labelname.split('Q')[1]).split('-')[0])
mtotal = 1.35*(1+q)
if 'Q2' not in labelname:continue
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
#Compute GW frequency
phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
#Limit the data to a small time interval during early inspiral say 200-350M
t_low = 200
t_high = 400
idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
time_fft = t22[idx_fft]
#Fit over this data
gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
plt.subplot(121)
plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Psi4 Frequency Oscillations')
plt.xlim(100,620)
#plt.ylim(-0.005,0.005)
plt.ylim(-0.015,-0.01)
#plt.subplot(222)
#plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
#plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('Psi4 Frequency Oscillations')
#plt.xlim(250,350)
##plt.ylim(-0.015,-0.01)
#plt.ylim(-0.1,0.1)
#Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
idx_fft = np.intersect1d(np.where(t22>t_low+10), np.where(t22<t_high-30))
time_fft = t22[idx_fft]
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
osc_fft = np.fft.fft(filtered)
fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)
plt.subplot(122)
plt.plot(fft_freq_hz, np.abs(osc_fft.real+1.j*osc_fft.imag),ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$FFT Freq$')
plt.ylabel('FFT(Psi4 Frequency)')
plt.axvline(x=1880)
plt.xlim(0,3000)
#plt.ylim(-3,2.5)
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
# Q3 GW frequency oscillations
from scipy import signal
plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
q = float((labelname.split('Q')[1]).split('-')[0])
mtotal = 1.35*(1+q)
if 'Q3' not in labelname:continue
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
#Compute GW frequency
phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
#Limit the data to a small time interval during early inspiral say 200-300M
t_low = 200
t_high = 400
idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
time_fft = t22[idx_fft]
#Fit over this data
gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
plt.subplot(121) #221)
plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Psi4 Frequency Oscillations')
plt.xlim(200,400)
#plt.ylim(-0.005,0.005)
plt.ylim(-0.015,-0.01)
#plt.subplot(222)
#plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
#plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('Psi4 Frequency Oscillations')
#plt.xlim(220,380)
##plt.ylim(-0.015,-0.01)
#plt.ylim(-0.1,0.1)
#Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
idx_fft = np.intersect1d(np.where(t22>t_low+30), np.where(t22<t_high-30))
time_fft = t22[idx_fft]
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
osc_fft = np.fft.fft(filtered)
fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)
plt.subplot(122) #223)
plt.plot(fft_freq_hz, np.abs(osc_fft.real),ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$FFT Freq$')
plt.ylabel('FFT(Psi4 Frequency)')
#plt.axvline(x=0.028)
plt.xlim(500,2000)
#plt.ylim(-3,2.5)
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
# Q5 GW frequency oscillations
from scipy import signal
plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
q = float((labelname.split('Q')[1]).split('-')[0])
mtotal = 1.35*(1+q)
if 'Q5' not in labelname:continue
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
#Compute GW frequency
phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
#Limit the data to a small time interval during early inspiral say 200-300M
t_low = 200
t_high = 400
idx_fft = np.intersect1d(np.where(t22>t_low), np.where(t22<t_high))
time_fft = t22[idx_fft]
#Fit over this data
gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
tfreq_fit = np.arange(t_low, t_high, t22[1]-t22[0])
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
plt.subplot(121) # 221)
plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Psi4 Frequency Oscillations')
plt.xlim(t_low, t_high)
#plt.ylim(-0.005,0.005)
plt.ylim(-0.015,-0.01)
#plt.subplot(222)
#plt.plot(time_fft, gw_osc,ls=ls[i], c=color[i], label=labelname)
#plt.plot(time_fft, filtered,ls='--', c='k', label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('Psi4 Frequency Oscillations')
#plt.axvline(x=481)
#plt.axvline(x=511)
#plt.xlim(t_low+20,t_high-20)
##plt.ylim(-0.015,-0.01)
#plt.ylim(-0.1,0.1)
#Remove the first and last 20M of the fit to avoid side effects of polynomial fitting
idx_fft = np.intersect1d(np.where(t22>t_low+30), np.where(t22<t_high-30))
time_fft = t22[idx_fft]
gw_osc = (gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft)
#Apply high pass filter to remove noise
sos = signal.butter(10, 6000*(mtotal/s1_in_Msun), 'lp', output='sos', )
filtered = signal.sosfilt(sos, gw_osc)
osc_fft = np.fft.fft(filtered)
fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
fft_freq_hz = fft_freq/(mtotal/s1_in_Msun)
plt.subplot(122) #223)
plt.plot(fft_freq_hz, np.abs(osc_fft.real),ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$FFT Freq$')
plt.ylabel('FFT(Psi4 Frequency)')
#plt.axvline(x=0.028)
#plt.xlim(0,0.2)
plt.xlim(500,2000)
#plt.ylim(-3,2.5)
plt.tight_layout()
plt.legend()
plt.show()
plt.close()
In [ ]:
#Strain vs Time
color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
j=0
fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
#filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'Strain/Strain_l2_m2.txt')
if not os.path.exists(filepath):continue
massratio = int((labelname.split('-')[0]).split('Q')[1])
if massratio==2:
j=1
if massratio==3:
j=2
if massratio==5:
j=3
t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
strain = np.sqrt(str_re22**2 + str_im22**2)
max_amp_idx = np.where(strain==np.amax(strain))
t22_ret = t22-R
plt.subplot(1,3,j)
plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
#plt.plot(t22_ret - t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2),c=color[i])
plt.xlabel(r'$(T - R_{ex}/M_{0}$')
plt.ylabel(r'$r \;h(t)\;/M_{0}$')
plt.xlim(-100, 100)
#plt.ylim(1e-28, 1e-23)
plt.legend()
plt.tight_layout()
plt.legend()
plt.grid(True)
fig.suptitle('Figure 4.1: Strain h(t) Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()
#Strain vs Time
color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = int((labelname.split('-')[0]).split('Q')[1])
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
#filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'Strain/Strain_l2_m2.txt')
if not os.path.exists(filepath):continue
t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
strain = np.sqrt(str_re22**2 + str_im22**2)
str_phase = np.unwrap(np.arctan2(str_im22, str_re22))
str_freq = np.gradient(str_phase)/np.gradient(t22)/np.pi/2.
max_amp_idx = np.where(strain==np.amax(strain))
M_in_s1 = (1+massratio)*1.35/s1_in_Msun
#print(labelname,t22[np.argmax(t22>150)],str_freq[np.argmax(t22>150)], str_freq[np.argmax(t22>150)]/M_in_s1)
str_physunits = (strain/0.1)*2.4*1e-22*(1+massratio)*1.35/5.
t22_src = t22
plt.subplot(121)
#plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i])
plt.plot(t22 - t22[max_amp_idx], strain,c=color[i], label=labelname)
#plt.plot((t22 - t22[max_amp_idx])*M_in_s1, str_physunits,c=color[i], label=labelname)
plt.xlabel(r'$(T - R_{ex}/M_{0}$')
plt.ylabel(r'$r \;h(t)\;/M_{0}$')
plt.xlim(-100, 100)
#plt.xlim(-0.01, 0.003)
#plt.ylim(1e-28, 1e-23)
plt.legend()
plt.grid(True)
plt.subplot(122)
#plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i])
plt.plot(t22_src - t22_src[max_amp_idx],str_freq ,c=color[i], label=labelname)
plt.xlabel(r'$(T - R_{ex})/M_{0}$')
plt.ylabel(r'$ Freq$')
plt.xlim(-100, 100)
plt.ylim(0,0.1)
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.legend()
fig.suptitle('Figure 4.2: Strain h(t) Amplitude Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()
In [ ]:
#QNM Fitting - Thesis Figure
#Psi4 plots: BBH vs NSBH
from collections import defaultdict
i=0
color = ['k','b', 'orangered','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
lm_modes = {'2,2':(2,2), '2,1':(2,1), '3,3':(3,3), '3,2':(3,2)}#, '4,4':(4,4), '4,3':(4,3)}
y_amp_max = [1e-4, 5e-4, 1e-4, ]
y_phase_max = [120, 150,50,250]
y_phase_min = [-50, 20,-50,50]
qnm_tau_dict = defaultdict(dict)
qnm_freq_dict = defaultdict(dict)
fig = plt.figure(figsize=(14,16))
#fig = plt.figure(figsize=(14,20))
nsbh_iter = 0
for lm_idx, lm_key in enumerate(sorted(lm_modes.keys())):
lm=lm_modes[lm_key]
#print("\nModes: l=%d, m=%d \n"%(lm[0], lm[1]))
#print("Simulation \t QNM Tau \t QNM Freqeuncy")
#print("-"*50)
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
nsbh_iter+=1
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l%d_m%d_r70.00.asc'%(lm[0],lm[1]))
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
phase = -1*np.unwrap(np.arctan2(psi4_im22,psi4_re22))
#Compute QNM Modes
#For QNM - A = constant*exp(-t\Ï„) and phase = omega*t+phi_0
#Find the fitting time interval
max_amp_idx = np.argmax(amp==np.amax(amp))
qnm_begin_idx = np.amin(np.where(t22 > t22[max_amp_idx]+40))
qnm_end_idx = np.amin(np.where(t22 > t22[max_amp_idx]+70))
#Create the fits
qnm_t_int = t22[qnm_begin_idx:qnm_end_idx]
qnm_amp_int = amp[qnm_begin_idx:qnm_end_idx]
qnm_phase_int = phase[qnm_begin_idx:qnm_end_idx]
logamp_fit_coeff = np.polyfit(qnm_t_int, np.log(qnm_amp_int),1)
logamp_fit = np.poly1d(logamp_fit_coeff)
phase_fit_coeff = np.polyfit(qnm_t_int, qnm_phase_int,1)
phase_fit = np.poly1d(phase_fit_coeff)
logamp_fit_slope = logamp_fit_coeff[0]
qnm_tau = -1./logamp_fit_slope
qnm_freq = phase_fit_coeff[0]/2./np.pi
qnm_tau_dict[lm_key][labelname] = qnm_tau
qnm_freq_dict[lm_key][labelname] = qnm_freq
#print("%10s \t %g \t %g "%(labelname, qnm_tau, qnm_freq))
#if lm[1]==1 or lm[1]==4: continue
#print(lm_idx, 2*lm_idx+1, 2*lm_idx+2)
plt.subplot(4,2,2*lm_idx+1)
amp_plot = plt.semilogy(t22-t22[max_amp_idx], amp,ls='-', c=color[i],label=labelname[:2])
#plt.plot(t22[max_amp_idx:]-t22[max_amp_idx], amp_fit(t22[max_amp_idx:]),ls='--', c=color[i])
#plt.plot(qnm_t_int-t22[max_amp_idx], logamp_fit(qnm_t_int),ls='-', c=color[i])
if lm_idx==2:plt.xlabel(r'$T/M$')
plt.ylabel(r'$|\psi_4|$')
plt.xlim(0,130)
#plt.ylim(1e-9,4e-4)
plt.legend()
plt.grid(True)
plt.subplot(4,2,2*lm_idx+2)
plt.plot(t22 - t22[max_amp_idx], phase,ls='-', c=color[i], label=labelname[:2])
#plt.plot(qnm_t_int - t22[max_amp_idx], phase_fit(qnm_t_int), ls='-', c=color[i])#, label=labelname+'Fit')
if lm_idx==2:plt.xlabel(r'$T/M$')
plt.ylabel(r'$\phi$')
plt.xlim(0,130)
plt.ylim(y_phase_min[lm_idx], y_phase_max[lm_idx])
plt.legend()
plt.grid(True)
#fig.suptitle('(%s) Mode'%lm_key, fontsize=16, y=1.08)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'QNMFitting.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
# Fit QNM Omega (from Berti's datafiles) using python fitting function
def fit_omega(a_bh, f1, f2, f3):
Momega_lmn = (f1 + f2*(1.-a_bh)**f3)
return Momega_lmn
def fit_Q(a_bh, q1, q2, q3):
Q_lmn = (q1 + q2*(1.-a_bh)**q3)
return Q_lmn
lmn = '210'
j, Momega_re, Momega_img = np.loadtxt('/nethome/bkhamesra3/Desktop/rsync/l2/n1l%sm%s.dat'%(lmn[0], lmn[1]), unpack=True, usecols=(0,1,2))
tau_by_M = -1./Mtau_inv
Q = Momega_re*tau_by_M/2.0
popt_omega, pcov_omega = curve_fit(fit_omega, j, Momega_re)
popt_tau, pcov_tau = curve_fit(fit_Q, j, Q)
a = 0.6
print(1./(Q[j==a]*2/Momega_re[j==a]), 1./(fit_Q(a, *popt_tau)*2/fit_omega(a, *popt_omega)))
plt.figure(figsize=(15,6))
plt.subplot(121)
plt.plot(j, Momega_re, label='Data')
plt.plot(j, fit_omega(j, *popt_omega), label='fit - Bhavesh')
plt.plot(j, fit_omega(j, f1_dict[lmn], f2_dict[lmn],f3_dict[lmn]), label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel(r'$\omega$')
plt.legend()
plt.subplot(122)
plt.plot(j, (Momega_re - fit_omega(j, *popt_omega))/Momega_re*100, label='fit - Bhavesh')
plt.plot(j,(Momega_re - fit_omega(j, f1_dict[lmn], f2_dict[lmn],f3_dict[lmn]))/Momega_re*100, label='fit - Berti')
plt.xlabel('a')
plt.ylabel(r'$\omega$ ')
plt.legend()
plt.show()
plt.close()
plt.figure(figsize=(15,6))
plt.subplot(121)
plt.plot(j, Q, label='Data')
plt.plot(j, fit_Q(j, *popt_tau), label='fit - Bhavesh')
plt.plot(j, fit_Q(j, q1_dict[lmn], q2_dict[lmn],q3_dict[lmn]), label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel('Q ')
plt.legend()
plt.subplot(122)
plt.plot(j, (Q - fit_Q(j, *popt_tau))/Q*100, label='fit - Bhavesh')
plt.plot(j,(Q - fit_Q(j, q1_dict[lmn], q2_dict[lmn],q3_dict[lmn]))/Q*100, label='fit -Berti Paper')
plt.xlabel('a')
plt.ylabel('Q % error')
plt.legend()
plt.show()
plt.close()
In [ ]:
#Strain - Mode Mixing - Thesis Figure
color = ['C0','C1','C2','r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
j=0
fig=plt.figure(figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
#filepath = os.path.join(datadir, 'Strain_l2_m2.txt')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
if 'Q5' not in labelname: continue
lm_modes = [(2,1), (2,2), (3,3), (4,4)]
for j,lm in enumerate(lm_modes):
filepath = os.path.join(datadir, 'Strain/Strain_l%d_m%d.txt'%(lm[0], lm[1]))
if not os.path.exists(filepath):continue
massratio = int((labelname.split('-')[0]).split('Q')[1])
if 'NS' in labelname:
plt.subplot(121)
ls='-'
else:
plt.subplot(122)
ls='--'
t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
strain = np.sqrt(str_re22**2 + str_im22**2)
max_amp_idx = np.where(strain==np.amax(strain))
t22_ret = t22-R
#plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
plt.semilogy(t22_ret -t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2), c=color[j],ls=ls,label='(%d, %d)'%(lm[0],lm[1]))
plt.xlabel(r'$(T - R_{ex})/M_{0}$')
plt.ylabel(r'$r \;h(t)\;/M_{0}$')
plt.xlim(t22_ret[0] -t22_ret[max_amp_idx]+150, 50)
plt.ylim(1e-3, 1e0)
plt.legend()
plt.grid(True)
plt.title(r'$%s - q=5$'%labelname[3:], y=1.01, fontsize=18)
plt.subplot(122)
lm_modes = [(2,1), (2,2), (3,3), (4,4)]
for j,lm in enumerate(lm_modes):
filepath = os.path.join(datadir, 'Strain/Strain_l%d_m%d.txt'%(lm[0], lm[1]))
if not os.path.exists(filepath):continue
t22, str_re22, str_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
strain = np.sqrt(str_re22**2 + str_im22**2)
max_amp_idx = np.where(strain==np.amax(strain))
t22_ret = t22-R
#plt.plot(t22_ret - t22_ret[max_amp_idx], str_re22, c=color[i], label=labelname)
plt.semilogy(t22_ret -t22_ret[max_amp_idx], np.sqrt(str_re22**2 + str_im22**2), ls='-', c=color[j])#,label='(%d, %d)'%(lm[0],lm[1]))
plt.xlabel(r'$(T - R_{ex})/M_{0}$')
plt.ylabel(r'$r \;h(t)\;/M_{0}$')
plt.xlim(t22_ret[0] -t22_ret[max_amp_idx]+150, 50)
plt.ylim(1e-3, 1e0)
plt.legend()
plt.grid(True)
#plt.title(r'$%s$'%labelname[3:], y=1.01)
plt.tight_layout()
plt.savefig(os.path.join(figdir,'Strain_ModeMixing.png'), dpi=100)
#fig.suptitle('Figure 4.1: Strain h(t) Comparison', fontsize=16, y=1.08)
plt.show()
plt.close()
Remnant BH Kicks¶
In [ ]:
#Remnant BH Kick
from scipy.constants import c as vc
plt.figure(figsize=(12,6))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
rad_arr = np.array((50,60,70,80,90,100,110,120,130,140))
Kmag_arr = []
for r in rad_arr:
filepath = os.path.join(datadir, 'psi4radAway_l6_r%d.00.asc'%r)
if not os.path.exists(filepath):continue
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
Kmag = (vc/1000)*(Kx**2 + Ky**2 + Kz**2)**0.5
Kmag_arr.append(Kmag[-1])
#print ('{:<16s} \t {:>1.2g} \t {:>1.2g} \t {:>1.2g} \t {:>1.4g}'.format(labelname,Kx[-1],Ky[-1],Kz[-1], Kmag[-1] ))
compare_nsbh[key]['Kick'] = Kmag_arr
compare_nsbh[key]['Det_radius'] = rad_arr
if len(Kmag_arr)<1: continue
plt.plot(compare_nsbh[key]['Det_radius'], compare_nsbh[key]['Kick'], label=labelname)
plt.xlabel('Radius/M')
plt.ylabel('Kick Velocity')
plt.legend()
plt.show()
plt.close()
Energy¶
In [ ]:
#Energy Plot:
i=0
plt.figure(figsize=(15,8))
color = [ 'k','b','orangered','k','b','orangered' ,'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls='--'
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls='-'
filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):
print('File note found for {}'.format(labelname))
t, Erad = np.loadtxt(filepath, unpack=True, usecols=(0,1))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Radiated Energy')
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Erad,ls=ls, c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Radiated Energy')
plt.xlim(-400,200)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
ls='-'
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls='--'
continue
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):continue
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Edot,ls=ls, c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Luminosity')
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Edot,ls=ls, c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Luminosity')
plt.xlim(-300,200)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
In [ ]:
#Radiated Energy - Thesis Figure
i=0
plt.figure(figsize=(15,8))
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')#'psi4radAway_l6_r70.00.asc')
if not os.path.exists(filepath):
print('File note found for {}'.format(labelname))
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Erad,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Radiated Energy')
plt.xlim(-400,200)
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Edot,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Luminosity')
plt.xlim(-150,100)
plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'RadiatedEnergy.png'), dpi=100)
plt.show()
plt.close()
Anuglar Momentum (Jz)¶
In [ ]:
#Angular Momentum Plots:
i=0
fig = plt.figure(figsize=(15,8))
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
ls = '-'
c = color[i-3]
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls='--'
c = color[i]
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
#filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
#t2, Erad2, Kx2, Ky2, Kz2, Edot2, Pxdot2, Pydot2, Pzdot2, Jxdot2, Jydot2, Jzdot, Jx2, Jy2, Jz2\
# = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
#plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
plt.xlabel(r'$T/M$')
plt.ylabel(r'$J_z$')
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
#plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
plt.xlabel(r'$T/M$')
plt.ylabel(r'$J_z$')
plt.xlim(-50,80)
plt.legend()
print('\n')
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Plot1.png'))
plt.show()
plt.close()
In [ ]:
#Angular Momentum Plots:
i=0
fig = plt.figure(figsize=(15,8))
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
ls = '-'
c = color[i-3]
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls='--'
c = color[i]
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
#filepath = os.path.join(datadir, 'psi4analysis_r70.00.asc')
#t2, Erad2, Kx2, Ky2, Kz2, Edot2, Pxdot2, Pydot2, Pzdot2, Jxdot2, Jydot2, Jzdot, Jx2, Jy2, Jz2\
# = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Jz,ls=ls, c=c, label=labelname)
#plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
plt.xlabel(r'$T/M$')
plt.ylabel(r'$J_z$')
plt.xlim(-150,300)
plt.legend()
plt.subplot(122)
if 'NSBH' in labelname:
plt.plot(t_bh-(t22[max_amp_idx]-70), sz_bh,ls=ls, c=c, label=labelname)
if 'BBH' in labelname:
#plt.plot(t_bh[t_bh>(t22[max_amp_idx]-70)]-(t22[max_amp_idx]-70), sz_bh[t_bh>(t22[max_amp_idx]-70)],ls=ls, c=c, label=labelname)
plt.axhline(y = sz_bh[-1], ls=ls, color=c, label=labelname)
#plt.plot(t2-t22[max_amp_idx], Jz2,ls=ls[i-1], c=color[i-1], label=labelname+'psi4analysis')
plt.xlabel(r'$(T - r)/M$')
plt.ylabel(r'$S_z$')
plt.xlim(-150,300)
#plt.ylim(0.4,0.5)
plt.legend()
#print(t[np.amin(np.where(t-(t22[max_amp_idx])>=300.0))], t_bh[np.amin(np.where(t_bh-(t22[max_amp_idx])>=300.0))])
print('%s \t \t %g \t %g \t %g'%(labelname, Jz[np.amin(np.where(t-(t22[max_amp_idx]-70)>=300.0))], \
sz_bh[np.amin(np.where(t_bh-(t22[max_amp_idx]-70)>=300.0))], \
Jz[np.amin(np.where(t-(t22[max_amp_idx])>=300.0))]+ sz_bh[np.amin(np.where(t_bh-(t22[max_amp_idx]-70)>=300.0))]))
print('\n')
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Plot2.png'))
plt.show()
plt.close()
In [ ]:
0.40391 - 0.406509
In [ ]:
#Angular Momentum Plots:
i=0
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
fig = plt.figure(figsize=(8,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls = '-'
c = color[i]
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_0.asc')
elif(compare_nsbh[key]['type']=='BBH'):
continue
datadir = compare_nsbh[key]['dir']
ihspin_filepath = os.path.join(datadir, 'ihspin_hn_3.asc')
ls='--'
c = color[i]
adm_filepath = os.path.join(datadir, 'adm_ejp_det_1.asc')
R_ex_adm = 100
t_adm, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm = np.loadtxt(adm_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))
#Read Psi4, Radiated Quantitiess and Spin data
psi4rad_filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
R_ex_gw = np.float((os.path.basename(psi4rad_filepath).split('.asc')[0]).split('_r')[1])
t_rad, Erad, Kx_rad, Ky_rad, Kz_rad, Edot_rad, Pxdot_rad, Pydot_rad, Pzdot_rad,\
Jxdot_rad, Jydot_rad, Jzdot_rad, Jx_rad, Jy_rad, Jz_rad\
= np.loadtxt(psi4rad_filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
t_bh, sx_bh, sy_bh, sz_bh = np.loadtxt(ihspin_filepath, unpack=True, usecols=(0,1,2,3))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
#Make spline fit of each and compute errors
Jz_adm_splinefit = splrep(t_adm[::2],Jz_adm[::2],s=0,k=3)
Jz_rad_splinefit = splrep(t_rad[::16],Jz_rad[::16],s=0,k=3)
sz_bh_splinefit = splrep(t_bh[::16],sz_bh[::16],s=0,k=3)
Jz_adm_fit = splev(t_adm,Jz_adm_splinefit, der=0)
Jz_rad_fit = splev(t_rad, Jz_rad_splinefit, der=0)
Sz_fit = splev(t_bh, sz_bh_splinefit, der=0)
t_start_fit = 70+R_ex_adm
Jz_adm_fiterr = (Jz_adm[t_adm>t_start_fit]-Jz_adm_fit[t_adm>t_start_fit])/Jz_adm[t_adm>t_start_fit]*100
t_Jzadm_err = t_adm[t_adm>t_start_fit]
t_maxerr_Jzadm = t_Jzadm_err[Jz_adm_fiterr==np.amax(Jz_adm_fiterr)]
t_minerr_Jzadm = t_Jzadm_err[Jz_adm_fiterr==np.amin(Jz_adm_fiterr)]
t_start_fit = 70+R_ex_gw
Jz_rad_fiterr = (Jz_rad[t_rad>t_start_fit]-Jz_rad_fit[t_rad>t_start_fit])/Jz_rad[t_rad>t_start_fit]*100
t_Jzrad_err = t_rad[t_rad>t_start_fit]
t_maxerr_Jzrad = t_Jzrad_err[Jz_rad_fiterr==np.amax(Jz_rad_fiterr)]
t_minerr_Jzrad = t_Jzrad_err[Jz_rad_fiterr==np.amin(Jz_rad_fiterr)]
#t_start_fit = 70
#Sz_fiterr = (sz_bh[t_bh>t_start_fit]-Sz_fit[t_bh>t_start_fit])/sz_bh[t_bh>t_start_fit]*100
#t_Sz_err = t_bh[t_bh>t_start_fit]
#t_maxerr_Sz = t_Sz_err[Sz_fiterr==np.amax(Sz_fiterr)]
#t_minerr_Sz = t_Sz_err[Sz_fiterr==np.amin(Sz_fiterr)]
print(labelname)
print("ADM Ang Mom Fitting Error: Max = %g %% at t = %g, Min = %g%% at t=%g"%(np.amax(Jz_adm_fiterr), t_maxerr_Jzadm,\
np.amin(Jz_adm_fiterr), t_minerr_Jzadm))
print("Rad Ang Mom Fitting Error: Max = %g %% at t = %g, Min = %g%% at t=%g"%(np.amax(Jz_fiterr), t_maxerr_Jz,\
np.amin(Jz_fiterr), t_minerr_Jz))
print("Spin Fitting Error:Max = %g %% at t = %g, Min = %g %% at t=%g"%(np.amax(Sz_fiterr), t_maxerr_Sz,
np.amin(Sz_fiterr), t_minerr_Sz))
#Make spline fit of each in COM frame (shift time to COM frame)
Jz_adm_splinefit = splrep(t_adm[::2]-R_ex_adm,Jz_adm[::2],s=0,k=3)
Jz_rad_splinefit = splrep(t_rad[::16]-R_ex_gw,Jz_rad[::16],s=0,k=3)
sz_bh_splinefit = splrep(t_bh[::16],sz_bh[::16],s=0,k=3)
#Compute fits for common time
Jz_adm_fit = splev(t_bh,Jz_adm_splinefit, der=0)
Jz_rad_fit = splev(t_bh, Jz_rad_splinefit, der=0)
Sz_fit = splev(t_bh, sz_bh_splinefit, der=0)
Jz_matter = Jz_adm_fit - (Jz_rad_fit + Sz_fit)
#t_common = np.intersect1d(np.around((t_bh + r_ex), decimals=1), np.around(t, decimals=1))
#idx_rad_gw = np.isin(np.around(t, decimals=1), t_common)
#idx_bh = np.isin(np.around(t_bh+r_ex, decimals=1), t_common)
#
#print(len(t_common), len(t_bh[idx_bh]) , len(t[idx_rad_gw]))
#if 'Q5' not in labelname: continue
#plt.subplot(121)
plt.plot(t_bh - (t22[max_amp_idx] -70), Jz_matter,ls='-', c=c, label=labelname)
#plt.plot(t_adm , Jz_adm_fit,ls='--', c=color[i-1], label=labelname+'Fitting')
plt.xlabel(r'$T /M $')
plt.ylabel(r'$J^{remining}_z$')
plt.ylim(0,1)
plt.xlim(-500,500)
plt.legend()
#plt.subplot(122)
#plt.plot(t_Jzadm_err , Jz_adm_fiterr,ls='-', c=c, label=labelname)
#
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel(r'$\Delta J_ADM(\%)$')
##plt.xlim(-50,100)
#plt.legend()
#fig.suptitle('Angular Momentum', y=1.08, fontsize=20)
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'AngMom_Jz_Matter_Plot3.png'))
plt.show()
plt.close()
In [ ]:
#BH Irreducible Mass Plots - Comparison of mass and spin of final black hole at 300M after merger
import matplotlib.gridspec as gridspec
color = ['k','b', 'orangered','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_irrmass_dict = {}
bbh_hznmass_dict = {}
print("System \t \t Init. Irr Mass Final Irr Mass \t NSBH - BBH \t Final Spin(a) \tFinal Hzn Mass NSBH - BBH \n ")
print('-'*115)
fig = plt.figure(constrained_layout=True, figsize=(16,8))
gs = gridspec.GridSpec(ncols=2, nrows=2, figure=fig)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])
ax3 = fig.add_subplot(gs[1,:])
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = str((labelname.split('Q')[1]).split('-')[0])
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 500.
bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
bbh_irrmass_bh1 = irrmass_ah1[0]
bbh_spinmag_bh1 = (sx1**2 + sy1**2 + sz1**2)**0.5
bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
#plt.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
#plt.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
#if 'Q2' in labelname:
# ax1.plot(t_ah1[:-100]-t_merge, irrmass_ah1[:-100], c=color[i], ls=ls[i],label=labelname)
# ax1.plot(t_ah3[530:]-t_merge, irrmass_ah3[530:], c=color[i], ls=ls[i])#,label=labelname)
#else:
# ax1.plot(t_ah1-t_merge, irrmass_ah1, c=color[i], ls=ls[i],label=labelname)
# ax1.plot(t_ah3-t_merge, irrmass_ah3, c=color[i], ls=ls[i])#,label=labelname)
ax1.axhline(y=irrmass_ah3[-1],c=color[i], ls='--')
#ax2.plot(t_ah1-t_merge, np.gradient(irrmass_ah1)/np.gradient(t_ah1), c=color[i], ls=ls[i],label=labelname)
#ax2.plot(t_ah3-t_merge, np.gradient(irrmass_ah3)/np.gradient(t_ah3), c=color[i], ls=ls[i])#,label=labelname)
#Set all the simulations to same final time before comparing mass and spins
area_ah3 = area_ah3 [t_ah3 < t_comp]
irrmass_ah3 = irrmass_ah3 [t_ah3 < t_comp]
arealrad_ah3 = arealrad_ah3[t_ah3 < t_comp]
t_ah3 = t_ah3 [t_ah3 < t_comp]
sx3 = sx3 [t_ih3 < t_comp]
sy3 = sy3 [t_ih3 < t_comp]
sz3 = sz3 [t_ih3 < t_comp]
t_ih3 = t_ih3[t_ih3 < t_comp]
bbh_irrmass_bh3 = irrmass_ah3[-1]
bbh_spinmag_bh3 = (sx3**2 + sy3**2 + sz3**2)**0.5
bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
bbh_irrmass_dict[massratio] = bbh_irrmass_bh3
bbh_hznmass_dict[massratio] = bbh_hznmass_bh3
print("%-10s \t %g \t %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
bbh_hznmass_bh3, bbh_hznmass_bh3 - bbh_hznmass_bh3))
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 500.
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
if 'Q5' in labelname:irrmass_splinefit = splrep(t_ah[::16],irrmass_ah[::16],s=0)
else :irrmass_splinefit = splrep(t_ah[::128],irrmass_ah[::128],s=0)
#Compute mass growth rate on every 8th iteration to reduce the high frequency noise in the derivatives
#Truncation and round off errors leads to high frequency noise which gets amplified on computing the derivatives
#Computing the derivative on a larger time step averages out these effect to zero)
mass_growth_rate_fit = splev(t_ah, irrmass_splinefit, der=1)
mass_growth_rate = np.gradient(irrmass_ah[::32])/np.gradient(t_ah[::32])
#plt.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
ax1.plot(t_ah-t_merge, irrmass_ah, c=color[i-3], ls='-',label=labelname[:2])
#ax1.plot(t_ah-t_merge, splev(t_ah, irrmass_splinefit), color[i], ls='--')
ax2.semilogy((t_ah-t_merge), mass_growth_rate_fit, c=color[i-3], ls='-',label=labelname[:2])
#ax2.semilogy((t_ah-t_merge)[::32], mass_growth_rate, c=color[i], ls='--',label=labelname[:2])
#Set all the simulations to same final time before comparing mass and spins
area_ah = area_ah [t_ah<t_comp]
irrmass_ah = irrmass_ah [t_ah<t_comp]
arealrad_ah = arealrad_ah[t_ah<t_comp]
t_ah = t_ah [t_ah<t_comp]
sx = sx [t_ih<t_comp]
sy = sy [t_ih<t_comp]
sz = sz [t_ih<t_comp]
t_ih = t_ih[t_ih<t_comp]
nsbh_irrmass_bh1 = irrmass_ah[0]
nsbh_spinmag_bh1 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
nsbh_irrmass_bh3 = irrmass_ah[-1]
nsbh_spinmag_bh3 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
print("%-1s \t %g \t %g \t\t %g \t %g \t %g \t %g"%(labelname, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
(nsbh_irrmass_bh3-bbh_irrmass_dict[massratio])/bbh_irrmass_dict[massratio], nsbh_dimspin_bh3,\
nsbh_hznmass_bh3, (nsbh_hznmass_bh3-bbh_hznmass_dict[massratio])/bbh_hznmass_dict[massratio]))
ax1.set_xlabel(r'$T/M$', fontsize=20)
ax1.set_ylabel(r'$ M_h$', fontsize=20)
ax2.set_xlabel(r'$T/M$', fontsize=20)
ax2.set_ylabel(r'$\dot{M}_h$', fontsize=20)
ax1.set_xlim(-200,300)
ax2.set_xlim(-50,100)
#ax1.set_ylim(-2,5)
ax2.set_ylim(2e-5,0.1)
ax1.legend(loc=4)
ax2.legend()
#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)
ax1.grid(True)
ax2.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
plt.show()
plt.close()
In [ ]:
#Rate of change of Angular Momentum
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
ls = '-'
c = color[i-3]
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
ls = '--'
c = color[i]
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Jzdot,ls=ls, c=c, label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel(r'$dJ_z/dt$')
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Jzdot,ls=ls, c=c, label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel(r'$dJ_z/dt$')
plt.xlim(-50,80)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
In [ ]:
#Remnant Kicks:
i=0
plt.figure(figsize=(15,8))
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Kx,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Kx')
plt.xlim(-50,80)
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Ky,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Ky')
plt.xlim(-50,80)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())): #[0:3:2])):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'psi4radAway_l6_r70.00.asc')
t, Erad, Kx, Ky, Kz, Edot, Pxdot, Pydot, Pzdot, Jxdot, Jydot, Jzdot, Jx, Jy, Jz\
= np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Pxdot,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel(r'$dP_x/dt$')
plt.xlim(-50,80)
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Pydot,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel(r'$dP_y/dt$')
plt.xlim(-50,80)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
ADM Quantities¶
In [ ]:
#ADM Energy Plot:
i=0
plt.figure(figsize=(15,8))
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
continue
datadir = compare_nsbh[key]['dir']
filepath = os.path.join(datadir, 'adm_ejp_det_3.asc')
t, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Eadm,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('ADM Energy')
plt.legend()
plt.subplot(122)
plt.plot(t-t22[max_amp_idx], Eadm,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('ADM Energy')
plt.xlim(-50,80)
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
In [ ]:
#ADM Angular momentum Plot:
i=0
plt.figure(figsize=(10,8))
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'adm_ejp_det_2.asc')
t, Eadm, Px_adm, Py_adm, Pz_adm, Jx_adm, Jy_adm, Jz_adm = np.loadtxt(filepath, unpack=True, usecols=(0,1,2,3,4,5,6,7))
psi4_filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(psi4_filepath):
psi4_filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
t22, psi4_re22, psi4_im22 = np.loadtxt(psi4_filepath, unpack=True, usecols=(0,1,2))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
#plt.subplot(231)
#plt.plot(t-t22[max_amp_idx], Jx_adm,ls=ls[i], c=color[i], label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('$J_x^{ADM}$')
#plt.legend()
#
#plt.subplot(232)
#plt.plot(t-t22[max_amp_idx], Jx_adm,ls=ls[i], c=color[i], label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('$J_x^{ADM}$')
#plt.xlim(-50,80)
#plt.legend()
#
#plt.subplot(233)
#plt.plot(t-t22[max_amp_idx], Jy_adm,ls=ls[i], c=color[i], label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('$J_y^{ADM}$')
#plt.legend()
#
#plt.subplot(234)
#plt.plot(t-t22[max_amp_idx], Jy_adm,ls=ls[i], c=color[i], label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('$J_y^{ADM}$')
#plt.xlim(-50,80)
#plt.legend()
#plt.subplot(121)
plt.plot(t-t22[max_amp_idx], Jz_adm,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$T/M$')
plt.ylabel('$J_z^{ADM}$')
plt.legend()
#plt.subplot(122)
#plt.plot(t-t22[max_amp_idx], Jz_adm,ls=ls[i], c=color[i], label=labelname)
#plt.xlabel(r'$Time - T_{max}$')
#plt.ylabel('$J_z^{ADM}$')
#plt.xlim(-500,500)
#plt.legend()
plt.tight_layout()
plt.savefig(os.path.join(figdir, 'ADM_AngMom_Jz.png'))
plt.show()
plt.close()
Density Oscillations¶
The star is initially not quite stable with strong oscillations in cenral density. The density changes by >25% compared to the initial value. The oscillations do not reduce much continuing all the way to the merger. However, from the above analysis, it appears that these oscillations do not appear to affect the dynamics and gravitational wave signature significantly to the first order. A more rigorous comparison study to test the effects of stellar oscillations is required.
IGR vs Whisky - In both the simulations, the oscillations in the stellar density are of the same order and frequency. However, near the merger, the disruption of star and process of formation of accretion disks have significant variations. But, once the disk starts to settle, the density profile looks quite similar (check the density snapshots) in both the codes with final densities ~$10^{13}-10^{14} g/cm^3$.
The complete density and velocity evolutions can be found here
- Whisky - https://github.gatech.edu/pages/bkhamesra3/bkhamesra3.github.io/research/nsbh/Separation_9M/NSBH_q2_D9_a0_C0.145_M105_Whisky/index.html
- IGR - https://github.gatech.edu/pages/bkhamesra3/bkhamesra3.github.io/research/nsbh/Separation_9M/NSBH_q3_D9_a0_C0.145_M70_IGR/index.html
SH_P1 - Density evolutions not provided in corrected paper for q=2
Some Tests to try -
- Compare Boosted TOV with Moving TOV obtained using BowenID
- Do a full numerical comparison with actual simulation dataset from another independent code such as LORENE.
In [ ]:
#Density Plots
i=0
filename = 'MinSearch0.asc'
fig, ax1 = plt.subplots(1,1,figsize=(15,8))
#left, bottom, width, height = [0.62, 0.57, 0.35, 0.4]
#ax2 = fig.add_axes([left, bottom, width, height])
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
print("Compactness \t Oscillation Amplitude \t Relative change (%)")
for i, direc in enumerate(sorted(compare_nsbh.keys())):
if compare_nsbh[direc]['type']=='BBH':
continue
labelname = compare_nsbh[direc]['label']
datadir = os.path.join(compare_nsbh[direc]['dir'], 'Summary/Data')
filepath = os.path.join(datadir, filename)
t, max_rho = np.loadtxt(filepath, unpack=True, usecols=(1,5))
idx_stablestar = np.intersect1d(np.where(t>100), np.where(t<120))
err = (np.amax(max_rho[idx_stablestar])-np.amin(max_rho[idx_stablestar]))/np.amin(max_rho[idx_stablestar])
if 'HR' not in labelname:
print("%s \t \t %g \t \t %g"%(labelname, np.amax(err*np.amin(max_rho[idx_stablestar])), np.amax(err)*100))
ax1.plot(t, max_rho/max_rho[0], ls=ls[i],c=color[i], label=labelname)
ax1.set_xlabel('Time')
ax1.set_ylabel(r'$\rho_{max}/\rho_{init}$')
ax1.legend()
ax1.set_ylim(0.8, 1.5)
ax1.set_xlim(200,400)
ax1.axvline(x=230)
ax1.axvline(x=270)
#ax2.plot(t, max_rho/max_rho[0], ls=ls[i],c=c[i])
#ax2.set_xlim(0, 100)
#ax2.set_ylim(0.9, 1.3)
#ax2.xaxis.tick_top()
#plt.xticks([])
#plt.yticks([])
i+=1
plt.tight_layout()
plt.ylim(0,2)
#plt.xlim(0,100)
#plt.savefig()
plt.show()
plt.close()
In [ ]:
#Combine density and velocity h5 files
import subprocess
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
else:
continue
vel_h5 = sorted(glob.glob(os.path.join(datadir, 'vel???.xy.h5')))
if len(vel_h5)<1:
velx_h5 = [os.path.join(datadir, 'velx.xy.h5')]
vely_h5 = [os.path.join(datadir, 'vely.xy.h5')]
velz_h5 = [os.path.join(datadir, 'velz.xy.h5')]
vel_h5 = velx_h5 + vely_h5 + velz_h5
hdf5_merge = ["/nethome/bkhamesra3/Cactus_branches/Cactus_2019/exe/ETK/hdf5_merge"]
rho_h5 = [os.path.join(datadir, 'rho.xy.h5')]
h5_files = rho_h5+vel_h5
process_string = hdf5_merge+h5_files+[os.path.join(datadir,'rho_vel.xy.h5')]
#print process_string
combine_h5 = subprocess.check_output(process_string) #, stdout=subprocess.PIPE)
In [ ]:
#Density snapshots with velocity profile
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
else:
continue
print("%s: Density - Veocity Profile Snapshots"%labelname)
mpl.rc('axes', labelsize=16, grid=False)
varname = 'density'
filepath = os.path.join(datadir, 'rho_vel.xy.h5')
zoom=1
varlim=None
annotate_grids=False
annotate_sphere=True
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
if os.path.exists(bhdiag):
it_bh,time_bh,x_bh,y_bh,z_bh, rad_bh = np.loadtxt(bhdiag, unpack=True, usecols=(0,1,2,3,4,7))
ds_allit = yt.load(os.path.join(datadir, 'rho_vel.xy.h5'))
iterations = ds_allit.get_all_iterations()
#Only look at merger snapshots
#iterations = iterations[int(3*len(iterations)/4):]
freq = int(len(iterations)/3)
if labelname=='NSBH-Whisky': init_idx = 2
else: init_idx=1
for j in iterations:
if os.path.exists(bhdiag):
time_j = time_bh[it_bh==j]
xj,yj,zj = x_bh[it_bh==j],y_bh[it_bh==j], z_bh[it_bh==j]
rad_j = rad_bh[it_bh==j]
ds = yt.load(filepath, iteration=j)
min_var, max_var = ds.all_data().quantities.extrema(varname)
var = yt.SlicePlot(ds, 'z', varname, axes_unit='code_length', origin='native', width=60, center = [0,0,0]) #center=[xnsj, ynsj, znsj])
if annotate_grids: var.annotate_grids()
var.set_cmap(varname, 'dusk') #'viridis')
var.annotate_velocity(plot_args={"headwidth": 2,'color':'white'})
if len(time_j)>0:var.annotate_text((0.45,0.9), 't = %.2f'%time_j, coord_system='figure', text_args={'color':'white'})
else:var.annotate_text((0.45,0.9), 'it = %d'%j, coord_system='figure', text_args={'color':'white'})
if os.path.exists(bhdiag) and annotate_sphere:
try:
var.annotate_sphere([xj,yj,zj], radius=(rad_j, 'code_length'), circle_args={'facecolor':'black','fill':0, 'edgecolor':'white'})
except YTPlotCallbackError:
print("Annotation Failed")
try:
var.zoom((zoom))
var.set_font_size(20)
var.set_xlabel('x (M)')
var.set_ylabel('y (M)')
if varname=='density':
var.set_zlim('density', 1.0e9, 1e16) #max_var)
elif varlim!=None:
var.set_zlim(varname, varlim[0], varlim[1])
else:
var.set_zlim(varname, 5e-2, max_var)
# This forces the ProjectionPlot to redraw itself on the AxesGrid axes.
plot = var.plots[varname]
cbar = plot.cb
# Finally, this actually redraws the plot.
var._setup_plots()
if varname=='density':
tick_labels = []
ticks=10**(np.arange(12,19))
for i in range(12,19):
tick_labels.append(r"$10^{%d}$"%(i))
# cbar.set_ticks(ticks.tolist())
# cbar.set_ticklabels(tick_labels)
# Finally, this actually redraws the plot.
#var._setup_plots()
except :
print("Annotation Failed")
# Ensure the colorbar limits match for all plots
var.show()
In [ ]:
#BH Horizon Mass Plots - Comparison of mass and spin of final black hole at 300M after merger
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_irrmass_dict = {}
bbh_hznmass_dict = {}
print("System \t \t Init. Irr Mass Final Irr Mass \t NSBH - BBH \t Final Spin(a) \tFinal Hzn Mass NSBH - BBH \n ")
print('-'*115)
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(15,6))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
massratio = str((labelname.split('Q')[1]).split('-')[0])
if (compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
bhdiag1 = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah1, area_ah1,irrmass_ah1, arealrad_ah1 = np.loadtxt(bhdiag1, unpack=True, usecols=(1,25,26,27))
ihspin1 = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih1, sx1, sy1, sz1 = np.loadtxt(ihspin1, unpack=True, usecols=(0,1,2,3))
bbh_irrmass_bh1 = irrmass_ah1[0]
bbh_spinmag_bh1 = (sx1**2 + sy1**2 + sz1**2)**0.5
bbh_hznmass_bh1 = (bbh_irrmass_bh1**2 + 0.25*bbh_spinmag_bh1[0]**2/bbh_irrmass_bh1**2)**0.5
bhdiag3 = os.path.join(datadir,'BH_diagnostics.ah3.gp')
t_ah3, area_ah3,irrmass_ah3, arealrad_ah3 = np.loadtxt(bhdiag3, unpack=True, usecols=(1,25,26,27))
ihspin3 = os.path.join(datadir,'ihspin_hn_2.asc')
t_ih3, sx3, sy3, sz3 = np.loadtxt(ihspin3, unpack=True, usecols=(0,1,2,3))
#plt.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
#plt.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
if 'Q2' in labelname:
ax1.plot(t_ah1[:-100]-t_merge, irrmass_ah1[:-100]-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
ax1.plot(t_ah3[530:]-t_merge, irrmass_ah3[530:]-irrmass_ah1[0], c=color[i], ls=ls[i])#,label=labelname)
else:
ax1.plot(t_ah1-t_merge, irrmass_ah1-irrmass_ah1[0], c=color[i], ls=ls[i],label=labelname)
ax1.plot(t_ah3-t_merge, irrmass_ah3-irrmass_ah1[0], c=color[i], ls=ls[i])#,label=labelname)
#ax2.plot(t_ah1-t_merge, np.gradient(irrmass_ah1)/np.gradient(t_ah1), c=color[i], ls=ls[i],label=labelname)
#ax2.plot(t_ah3-t_merge, np.gradient(irrmass_ah3)/np.gradient(t_ah3), c=color[i], ls=ls[i])#,label=labelname)
#Set all the simulations to same final time before comparing mass and spins
area_ah3 = area_ah3 [t_ah3 < t_comp]
irrmass_ah3 = irrmass_ah3 [t_ah3 < t_comp]
arealrad_ah3 = arealrad_ah3[t_ah3 < t_comp]
t_ah3 = t_ah3 [t_ah3 < t_comp]
sx3 = sx3 [t_ih3 < t_comp]
sy3 = sy3 [t_ih3 < t_comp]
sz3 = sz3 [t_ih3 < t_comp]
t_ih3 = t_ih3[t_ih3 < t_comp]
bbh_irrmass_bh3 = irrmass_ah3[-1]
bbh_spinmag_bh3 = (sx3**2 + sy3**2 + sz3**2)**0.5
bbh_hznmass_bh3 = (bbh_irrmass_bh3**2 + 0.25*bbh_spinmag_bh3[-1]**2/bbh_irrmass_bh3**2)**0.5
bbh_dimspin_bh3 = bbh_spinmag_bh3[-1]/bbh_hznmass_bh3**2
bbh_irrmass_dict[massratio] = bbh_irrmass_bh3
bbh_hznmass_dict[massratio] = bbh_hznmass_bh3
print("%-10s \t %g \t %g \t\t %g \t\t %g \t %g \t %g "%(labelname,bbh_irrmass_bh1, bbh_irrmass_bh3, \
bbh_irrmass_bh3 - bbh_irrmass_bh3, bbh_dimspin_bh3,\
bbh_hznmass_bh3, bbh_hznmass_bh3 - bbh_hznmass_bh3))
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
t_comp = t_merge + 300.
#if 'Q2' in labelname:
# datadir = os.path.join('/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_2/NSBH_q2_D9_a0_C0.145_M105_QLM/Summary/Data/')
bhdiag = os.path.join(datadir,'BH_diagnostics.ah1.gp')
t_ah, area_ah,irrmass_ah, arealrad_ah = np.loadtxt(bhdiag, unpack=True, usecols=(1,25,26,27))
ihspin = os.path.join(datadir,'ihspin_hn_0.asc')
t_ih, sx, sy, sz = np.loadtxt(ihspin, unpack=True, usecols=(0,1,2,3))
mass_growth_rate = np.gradient(irrmass_ah[::8])/np.gradient(t_ah[::8])
#plt.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
ax1.plot(t_ah-t_merge, irrmass_ah-irrmass_ah[0], c=color[i], ls=ls[i],label=labelname)
ax2.plot((t_ah-t_merge)[::8], mass_growth_rate, c=color[i], ls=ls[i],label=labelname)
#Set all the simulations to same final time before comparing mass and spins
area_ah = area_ah [t_ah<t_comp]
irrmass_ah = irrmass_ah [t_ah<t_comp]
arealrad_ah = arealrad_ah[t_ah<t_comp]
t_ah = t_ah [t_ah<t_comp]
sx = sx [t_ih<t_comp]
sy = sy [t_ih<t_comp]
sz = sz [t_ih<t_comp]
t_ih = t_ih[t_ih<t_comp]
nsbh_irrmass_bh1 = irrmass_ah[0]
nsbh_spinmag_bh1 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh1 = (nsbh_irrmass_bh1**2 + 0.25*nsbh_spinmag_bh1[0]**2/nsbh_irrmass_bh1**2)**0.5
nsbh_irrmass_bh3 = irrmass_ah[-1]
nsbh_spinmag_bh3 = (sx**2 + sy**2 + sz**2)**0.5
nsbh_hznmass_bh3 = (nsbh_irrmass_bh3**2 + 0.25*nsbh_spinmag_bh3[-1]**2/nsbh_irrmass_bh3**2)**0.5
nsbh_dimspin_bh3 = nsbh_spinmag_bh3[-1]/nsbh_hznmass_bh3**2
print("%-1s \t %g \t %g \t\t %g \t %g \t %g \t %g"%(labelname, nsbh_irrmass_bh1,nsbh_irrmass_bh3, \
(nsbh_irrmass_bh3-bbh_irrmass_dict[massratio])/bbh_irrmass_dict[massratio], nsbh_dimspin_bh3,\
nsbh_hznmass_bh3, (nsbh_hznmass_bh3-bbh_hznmass_dict[massratio])/bbh_hznmass_dict[massratio]))
ax1.set_xlabel(r'$T/M$')
ax1.set_ylabel(r'$\delta M_{h}$')
ax2.set_xlabel(r'$T/M$')
ax2.set_ylabel(r'$dM_{h}/dt$')
ax1.set_xlim(-200,300)
ax2.set_xlim(-50,100)
#ax1.set_ylim(-2,5)
#ax2.set_ylim(-0.5,1)
ax1.legend()
ax2.legend()
#plt.axhline(y=irrmass_ah[0], ls='--', c='k')
#plt.ylim(0.65, 0.92)
ax1.grid(True)
ax2.grid(True)
plt.tight_layout()
#plt.savefig(os.path.join(figdir, 'BH_MassPlots.png'), dpi=100)
#plt.show()
plt.close()
Neutron Star Masses¶
The mass of neutron star rest mass remains roughly constant till the disruption at which point it drastically decreases. After this point, since the matter is consumed by the black hole, the code cannot find the surface reliably and probably leads to unreliable results.
In [ ]:
#NS Mass Plots -
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
volint = os.path.join(datadir,'volume_integrals-GRMHD.asc')
if not(os.path.exists(volint)): continue
t, restmass = np.loadtxt(volint, unpack=True, usecols=(0,4))
if compare_nsbh[key]['symmetry']:
restmass= 2*restmass #accounting for symmetries
print("%s:Initial Mass = %g, Final mass = %g"%(labelname[:4], restmass[0],restmass[-1] ))
plt.plot(t-t_merge, restmass, c=color[i], ls=ls[i],label=labelname)
plt.xlabel(r'$(T - R_{ex})/M$')
plt.ylabel('Mass')
#plt.ylim(0., 2)
plt.grid(True)
plt.legend()
plt.show()
plt.close()
#NS Mass Plots -
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
plt.figure(figsize=(15,8))
k=0
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if(compare_nsbh[key]['type']=='NSBH'):
k+=1
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
t_merge = compare_nsbh[key]['t_merger_src']
volint = os.path.join(datadir,'volume_integrals-GRMHD.asc')
if not(os.path.exists(volint)): continue
t, restmass = np.loadtxt(volint, unpack=True, usecols=(0,4))
if compare_nsbh[key]['symmetry']:
restmass= 2*restmass #accounting for symmetries
print("%s:Initial Mass = %g, Final mass = %g"%(labelname[:4], restmass[0],restmass[-1] ))
plt.subplot(2,2,k)
plt.plot(t-t_merge, 100*(restmass[0] - restmass)/restmass[0], c=color[i], ls=ls[i],label=labelname)
plt.xlabel(r'$(T - R_{ex})/M$')
plt.ylabel(r'$\delta Mass/M_0(\%)$')
plt.ylim(0., 2)
plt.xlim(-t_merge, 50)
plt.grid(True)
plt.legend()
plt.show()
plt.close()
In [ ]:
#Accretion disk and dynamical Ejecta Rest Mass-
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_ah3_mass = 0
#subplt_num=1
print("System \t\t\t M(R>R_{AH})/M_* \t Accretion Disk Mass \t M_disk/M_* \t Dynamical Ejecta Mass/M_* \n ")
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
t_merge = compare_nsbh[key]['t_merger_src']
if (compare_nsbh[key]['type']=='BBH'):
continue
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
time, m_ns, m_bh3, m_0d7, m_1d4, m_20, m_35 = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5,6,7,8,9,10))
if compare_nsbh[key]['symmetry']:
m_ns = 2.*m_ns #accounting for symmetries
m_bh3 = 2.*m_bh3 #accounting for symmetries
m_0d7 = 2.*m_0d7 #accounting for symmetries
m_1d4 = 2.*m_1d4 #accounting for symmetries
m_20 = 2.*m_20 #accounting for symmetries
m_35 = 2.*m_35 #accounting for symmetries
if 'Q5' in labelname:tcap = 180
else:tcap=500
time = time-t_merge
idx_cap500 = np.where(time<=tcap)
m_ns = m_ns [idx_cap500]
m_bh3 = m_bh3 [idx_cap500]
m_0d7 = m_0d7 [idx_cap500]
m_1d4 = m_1d4 [idx_cap500]
m_20 = m_20 [idx_cap500]
m_35 = m_35 [idx_cap500]
time = time [idx_cap500]
print("%-15s \t %g \t\t %g \t\t %g\t \t %g "%(labelname, m_1d4[-1]/m_ns[0], m_1d4[-1]-m_35[-1], \
m_1d4[-1]/m_ns[0]-m_35[-1]/m_ns[0], m_35[-1]/m_ns[0] ))
plt.figure(figsize=(15,6))
plt.subplot(1, 2, 1)
#plt.plot(time, m_0d7, c='C0', ls=':', label=labelname+'R=0.7')
plt.plot(time, m_1d4/m_ns[0], c='C1', ls='-', label=labelname+' (R=1.5)')
plt.plot(time, m_35/m_ns[0], c='C2', ls='-.', label=labelname+' (R=35)')
plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
plt.xlabel('Time')
plt.ylabel(r'$M_{r>R}/M_*$')
plt.grid(True)
#plt.ylim(0.0, 0.3)
#plt.xlim(t_merge-100,1100) #time[-1]+10)
plt.legend()
plt.subplot(122)
#plt.semilogy(time, m_0d7, c='C0', ls=':', label=labelname+'R=0.7')
plt.semilogy(time, m_1d4/m_ns[0]-m_35[-1]/m_ns[0], c='C1', ls='-', label=labelname+' (R=1.5)')
plt.semilogy(time, m_35/m_ns[0], c='C2', ls='-.', label=labelname+' (R=35)')
plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
plt.xlabel('Time')
plt.ylabel(r'$M_{r>r_{R}}/M_*$')
#plt.ylim(0.9, 0.94)
#plt.xlim(t_merge-100,1100)
plt.legend()
plt.grid(True)
plt.show()
plt.close()
#Accretion disk and dynamical Ejecta Rest Mass-
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_ah3_mass = 0
In [ ]:
#Total Rest Mass outside BH Plots -
color = ['r', 'b', 'k','m','coral','forestgreen']
ls = ['-', '--', ':','-.','--',':']
bbh_ah3_mass = 0
#print("System \t\t\t Initial BH Mass \t Final BH Mass \t NSBH - BBH \n ")
#subplt_num=1
for i, key in enumerate(sorted(compare_nsbh.keys())):
#print("\n\n")
labelname = compare_nsbh[key]['label']
t_merge = compare_nsbh[key]['t_merger_src']
if (compare_nsbh[key]['type']=='BBH'):
continue
elif(compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
vi_grmhd = os.path.join(datadir,'volume_integrals-GRMHD.asc')
time, m_ns = np.loadtxt(vi_grmhd, unpack=True, usecols=(0,5))
if compare_nsbh[key]['symmetry']:
m_ns = 2.*m_ns #accounting for symmetries
totalrestmass = os.path.join(datadir,'TotalRestMass.asc')
if 'Q2' in labelname: totalrestmass = '/localdata2/bkhamesra3/simulations/Hive/NSBH/MassRatio_2/NSBH_q2_D9_a0_C0.145_M105_QLM/Summary/Data/TotalRestMass.asc'
data = np.loadtxt(totalrestmass, unpack=True, usecols=(0,1,2,4,6,8,10))
it = data[0]
print(it[2]-it[1])
if it[2]-it[1]<30:
data = (data.T[::128]).T
it=data[0]
print(it[2]-it[1])
it, time, m_0d7, m_1d4, m_15, m_25,m_35 = data
#print("%-15s \t %g \t \t %g \t %g"%(labelname, irrmass_ah[0],irrmass_ah[-1], irrmass_ah[-1]-bbh_ah3_mass ))
#plt.subplot(2, 1, subplt_num)
plt.figure(figsize=(15,8))
#plt.plot(time, m_0d7, c='C0', ls='solid',label=labelname+'(R=0.7)')
#plt.plot(time, m_1d4, c='C1', ls='--',label=labelname+'(R=1.4)')
#plt.plot(time, (m_6-m_1d4)/m_ns[0], c='C2', ls=':', label='M(R=6) - M(R=1.4)')
plt.plot(time, (m_25-m_1d4), c='C3', ls='-', label='M(R=25)- M(R=1.4)')
plt.plot(time, (m_35-m_25), c='C4', ls='-.', label='M(R=35)- M(R=25)')
plt.axvline(x=t_merge, ls='-', c='k', label='Time of merger')
#plt.yscale('log')
plt.xlabel('Time')
plt.ylabel(r'Rest Mass /$M_*$')
plt.title(labelname, fontsize=16)
#subplt_num+=1
plt.ylim(1e-6,1.2)
plt.xlim(0,1100)
plt.legend()
plt.show()
plt.close()
Rough Work¶
In [ ]:
#BBH Old and New comparisons
ylm1 = '/localdata2/bkhamesra3/simulations/BBH/Test/D9_q3.0_a0.0_m160/Ylm_WEYLSCAL4::Psi4r_l2_m2_r100.00.asc'
t, re, im = np.loadtxt(ylm1, unpack=True, usecols=(0,1,2))
amp = (re**2 + im**2)**0.5
plt.figure(figsize=(15,6))
plt.plot(t, re, c='C0', label='Old')
plt.plot(t, amp, c='C0',ls='--')
ylm1 = '/localdata2/bkhamesra3/simulations/BBH/Hive/BBH_q3_D9_a0_M250/BBH_q3_D9_a0_M250-all/Ylm_WEYLSCAL4::Psi4r_l2_m2_r100.00.asc'
t, re, im = np.loadtxt(ylm1, unpack=True, usecols=(0,1,2))
amp = (re**2 + im**2)**0.5
plt.plot(t, re, c='C1', label='New')
plt.plot(t, amp, c='C1',ls='--')
plt.xlim(700,950)
plt.legend()
plt.show()
plt.close()
In [ ]:
os.chdir('/localdata/bkhamesra3/research_localdata/UsefulScripts/Python_Notebooks/BowenID_Tests/NSBH/Simulation_FinalResults/MassRatio_3/')
In [ ]:
! jupyter nbconvert --to html NSBH_q3_D9_a0_C0.145_Comparison.ipynb
In [ ]:
! mv NSBH_q3_D9_a0_C0.145_Comparison.html /localdata2/bkhamesra3/website/github_pages/bkhamesra3.github.io/research/nsbh/Separation_9M/Comparisons/MassRatio3/
In [ ]:
a = np.array((1,2,3))
print((np.where(a>3) and a[-1]))
In [ ]:
# Q5 GW frequency oscillations
plt.figure(figsize=(15,8))
for i, key in enumerate(sorted(compare_nsbh.keys())):
labelname = compare_nsbh[key]['label']
if (compare_nsbh[key]['type']=='NSBH'):
datadir = os.path.join(compare_nsbh[key]['dir'], 'Summary/Data')
tmerge_nsbh = compare_nsbh[key]['t_merger_src']
elif(compare_nsbh[key]['type']=='BBH'):
datadir = compare_nsbh[key]['dir']
continue
filepath = os.path.join(datadir, 'Ylm_WEYLSCAL4::Psi4r_l2_m2_r70.00.asc')
if not os.path.exists(filepath):
filepath = os.path.join(datadir, 'mp_Psi4r_l2_m2_r70.00.asc')
q = float((labelname.split('Q')[1]).split('-')[0])
mtotal = 1.35*(1+q)
if 'Q5' not in labelname:continue
t22, psi4_re22, psi4_im22 = np.loadtxt(filepath, unpack=True, usecols=(0,1,2))
#Compute GW frequency
phi = np.unwrap(np.arctan2(psi4_im22,psi4_re22))
amp = np.sqrt(psi4_re22**2 + psi4_im22**2)
max_amp_idx = np.argmax(amp==np.amax(amp))
gw_freq = np.gradient(phi)/np.gradient(t22)/2./np.pi
#Limit the data to a small time interval during inspiral say 200-300M
idx_fft = np.intersect1d(np.where(t22>320), np.where(t22<420))
time_fft = t22[idx_fft]
gw_freq_fit = np.poly1d(np.polyfit(t22[idx_fft], gw_freq[idx_fft], deg=2))
tfreq_fit = np.arange(300, 500, t22[1]-t22[0])
osc_fft = np.fft.fft((gw_freq[idx_fft]-gw_freq_fit(time_fft))/gw_freq_fit(time_fft))
fft_freq = np.fft.fftfreq(np.size(time_fft), d=time_fft[2]-time_fft[1])
plt.subplot(121)
#plt.plot(t22, (gw_freq-gw_freq_fit(t22))/gw_freq_fit(t22),ls=ls[i], c=color[i], label=labelname)
plt.plot(t22, gw_freq,ls='-', c=color[i], label=labelname)
plt.plot(tfreq_fit, gw_freq_fit(tfreq_fit),ls='--', c='k', label=labelname)
plt.xlabel(r'$Time - T_{max}$')
plt.ylabel('Psi4 Frequency Oscillations')
plt.xlim(200,700)
plt.ylim(-0.02,-0.003)
#plt.ylim(-0.1,0.1)
plt.subplot(122)
plt.plot(fft_freq/(mtotal/s1_in_Msun), osc_fft.real,ls=ls[i], c=color[i], label=labelname)
plt.xlabel(r'$FFT Freq$')
plt.ylabel('FFT(Psi4 Frequency)')
plt.axvline(x=0.028)
plt.xlim(0,2000)
#plt.ylim(-3,2.5)
plt.legend()
plt.show()
plt.close()
In [ ]:
import h5py
h5file = '/localdata2/bkhamesra3/simulations/InjectionStudies/NSBH/Completed/NonSpinning/GT0014_NSBH.h5'
wf_data = h5py.File(h5file, 'r')
#$for atr_keys in wf_data.attrs.keys():
flow = wf_data.attrs['f_lower_at_1MSUN']
print(flow/8.1)
wf_data.close()